双层优化入门(4)—基于对偶变换的双层优化求解

        之前的博客介绍了双层优化的基本原理、以及如何使用KKT条件和智能优化算法求解双层优化问题，这篇博客将继续介绍如何通过对偶变换求解双层优化问题。

1.线性规划的对偶问题

        参考资料：

运筹学修炼日记：如何优雅地写出大规模线性规划的对偶_刘兴禄的博客

        首先说明一下，能应用强对偶定理求解的双层优化问题一定是线性规划问题，不能含有非线性目标函数和非线性约束，且不能含有0-1变量或者整数变量，如果有的话则需要考虑使用其他方法进行求解。

        用上面博客中提到的一个简单的线性规划来说明如何写出线性规划问题的对偶问题。

        原问题：

        对偶问题： 

        其中原问题为max问题，决策变量为x1和x2，对偶问题为min问题，决策变量为y1和y2，总结出的规律如下： 

        我们挨个来对比一下：

        1)原问题为max问题，对偶问题为min问题。

        2)原问题中有2个变量，因此对偶问题中有2个约束条件。

        3)原问题中变量都≥0，因此对偶问题中约束条件的符号也是≥0。

        4)原问题中目标函数的系数为2和3，因此对偶问题的约束条件中≥号右边的常数也为2和3。

        5)原问题中有2个约束条件，因此对偶问题有2个决策变量。

        6)原问题中约束条件的符号为≤，因此对偶问题中决策变量的取值都≥0。

        7)原问题中第1个约束条件中x1和x2的系数分别为5和4，因此对偶问题的约束条件中，变量y1在两个不同的约束条件中的系数分别为5和4。

        8)原问题中第2个约束条件中x1和x2的系数分别为2和3，因此对偶问题的约束条件中，变量y2在两个不同的约束条件中的系数分别为2和3。

        对于上面这种简单的优化问题，我们可以很容易写出其对偶问题，但如果原优化问题比较复杂，或者规模很大，就不太容易写出其对偶问题，即使写出来也容易出错。上面的博客中提到了用excel表格来进行对偶变换，这是一种比较好的方法。文献中另一种常见的做法是将优化问题先改写成紧凑的矩阵形式，然后使用矩阵变换就可以方便地表达原问题和对偶问题。

2.利用对偶变换求解双层优化

        通过上面的分析可知，经过对偶变换之后，优化问题的目标函数将从min变为max或者从max变为min，那么对于min-max或者max-min类型的双层优化问题，首先可以使用对偶变换将下层优化转为和上层优化方向一致，然后就可以将双层优化问题转为单层优化。

        再次总结双层优化可以使用对偶变换求解的条件：

        1.下层优化为线性规划，没有非线性项及0-1变量，整数变量。

        2.上下层优化的目标方向不一致，可以为min-max类型，也可以为max-min类型。

        3.上下层优化目标相同或者有一部分相同，这样下层优化经过对偶变换后就可以和上层优化进行合并。

        下面是一个简单线性双层优化的例子：

        首先这个双层优化问题只有两个变量，可以简单推断一下，这个双层优化问题的最优解中，x将取得其定义域内的最大值，y将取得其定义域内的最小值。当然猜想归猜想，还是要实际求解才知道是否正确。

        下面使用对偶变换求出该问题的最优解，验证我们的想法是否正确。我们可以把优化问题的目标函数进行拆分，然后把下层优化问题单独列出来：

        我们分析一下，下层优化问题共1个变量，5个约束条件，因此其对偶问题共有1个约束条件，5个变量。由于优化问题比较简单，我们直接写出其对偶问题： 

        进一步将这个对偶问题和上层优化问题合并可以得到： 

        其中Ωx表示原问题中x的可行域。 

        这样就把一个线性双层优化问题转换为了一个单层二次规划问题，可以使用Yalmip调用Cplex，Gurobi等求解器进行求解，代码如下：

%% 清空
clc
clear
close all
warning off
 
%% 目标函数和约束条件
x = sdpvar(1);
z = sdpvar(5,1);
y = sdpvar(1);
Constraints = [[-3 -1 -1 -1 1]*z <= 2 , z >= 0 , x >= 0 , y >= 0];
xy_Constraints = [-3*x+y <= -3 , 3*x+y <= 30 , 2*x-3*y >= -12 , x+y <= 14 , x >= 0 , y >= 0];
figure
plot(xy_Constraints,[],[],[],sdpsettings('plot.shade',0.1));
Constraints = [Constraints , boundingbox(xy_Constraints , [] , [x,y])];
objective = x + [-12-2*x x-14 3-3*x 3*x-30 0]*z;
ops = sdpsettings('verbose', 0 , 'solver', 'gurobi');
result = optimize(Constraints , -objective , ops);
 
%% 输出结果
disp(['最优解：x=',num2str(value(x)),',y=',num2str(value([-12-2*x x-14 3-3*x 3*x-30 0]*z/2))])
disp(['最优函数值=',num2str(value(objective))])

        运行结果：

      这个是变量x和y的取值范围。因为对偶变换之后，约束条件中都与变量z相关的项，没有对x进行约束，所以还需要从原来的上层优化问题中提取出变量x的取值范围，也就是上面提到的Ωx，Yalmip中可以直接使用boundingbox函数进行提取。

        由结果可知，当x取10，y取0时，整个双层优化有最优解10，这和我们根据常理进行的推断相同，x取得定义域内的最大值，y取得定义域内的最小值。