之前的博客：鲁棒优化入门（二）——基于matlab+yalmip求解鲁棒优化问题

去年发布了使用Yalmip工具箱求解鲁棒优化问题的博客之后，陆陆续续有朋友问我相关的问题，有人形容从学习这篇博客到求解论文中的鲁棒优化问题，就好像刚学会求导公式，就要去做高考压轴题，根本无从下手。为了解决这个问题，这篇博客将手把手地教会大家如何使用Matlab+ yalmip+cplex(当然其他的求解器比如gurobi也是可以的)求解论文中的鲁棒优化问题。为了具有拓展性，本文将选取两篇不同的参考文献，一个是标准的单阶段鲁棒优化问题，另一个是两阶段鲁棒优化问题，分别说明如何对鲁棒优化的文献进行分析，并借助Yalmip工具箱实现文献中单阶段和两阶段鲁棒优化的编程求解，这篇博客是单阶段鲁棒优化的实现。

1.单阶段鲁棒优化基本形式

如参考文献[1]中所述，标准的单阶段鲁棒优化的形式为：

其中x表示系统的状态(也可以认为是参数)，u表示人为控制变量(也就是决策变量)，w表示大自然的决策变量(也就是不确定变量)。万变不离其宗，几乎所有的单阶段鲁棒优化问题都可以写成这种形式，并采用相同的套路进行求解。因此，对于任意文献中的单阶段鲁棒优化问题，可以转为和上式统一的形式便于求解，具体转换思路如下：

1）确定优化问题中的决策变量有哪些，将其与变量u对应。

2）确定优化问题中的不确定变量是什么，将其与变量w对应。

3）确定优化问题中不确定集合的形式，并考虑是否可以直接使用Yalmip中的鲁棒优化模块进行求解。

4）确定优化问题的目标函数，将其与J(u,w)对应。

5）确定优化问题中的约束条件，是否包含非线性约束需要线性化。

2.选股优化问题Matlab+Yalmip编程实例

以一个经典的选股优化问题为例，说明如何使用Matlab实现编程求解单阶段鲁棒优化问题。

一共有n=150支股票，其中第i支股票的期望收益为pi，标准差为σi，用zi表示第i支股票和期望收益偏差的程度，xi表示第i只股票选取的比例，则这个问题可以描述为鲁棒优化问题：

2.1 鲁棒优化模块求解

按上面的思路逐步进行分析：

1）优化问题中的决策变量有哪些？

决策变量为第i只股票选取的比例xi。

2）优化问题中的不确定变量是什么？

不确定变量是第i支股票和期望收益偏差的程度zi。

3）确定优化问题中不确定集合的形式，并考虑是否可以直接使用Yalmip中的uncertain函数进行求解。

优化问题中不确定集合为多面体的形式，可以直接使用Yalmip的鲁棒优化功能进行求解。

4）确定优化问题的目标函数。

内层的目标函数为最小化股票的收益，即通过不确定变量zi的值找到最恶劣的场景；外层的目标函数为最大化股票收益，即通过确定决策变量xi的取值找到最恶劣场景下使股票收益最大化的方案。

5）确定优化问题中的约束条件。

约束条件包括各个股票选取的比例总和为1，各个股票选择的比例xi大于等于1以及不确定变量的约束等，不包含非线性约束，无需转换可以直接求解。

按上面的思路，Matlab编程如下：

%% 选股问题的鲁棒优化模块求解
clc
clear
close all
yalmip('clear')

%% 1.相关的参数
n  = 150;                                               % 股票的数量 
p  = 1.15 + 0.05/150*(1:n)';                            % 期望的收益
sigma = 0.05/450*sqrt(2*n*(n + 1)*(1:n)');              % 收益的偏差
gamma = 5;                                              % 不确定预算

%% 2.定义决策变量x和不确定变量z
z = sdpvar(n,1);                                        % 不确定变量z
x = sdpvar(n,1);                                        % 决策变量x

%% 3.确定目标函数
objective = -(p - sigma.*z)'*x;                         % 目标函数

%% 4.用uncertain函数构造不确定变量，并创建不确定集合
U = [uncertain(z) , z >= 0 , z <= 1 ,...                % 不确定集
    norm(z,Inf) <= 1 , norm(z,1) <= gamma];

%% 5.其他约束条件
Constraints = [sum(x) == 1 , x >= 0];                   % 其他约束条件

%% 6.求解优化问题
ops = sdpsettings('verbose', 3, 'solver', 'cplex' , 'showprogress' , 1);
sol = optimize(Constraints + U , objective , ops);      % 求解模型

%% 7.优化结果
x = value(x);                                           % 决策变量x取值
z = value(z);                                           % 不确定变量z的取值
plot(x)                                                 % 画出图像

运行结果如下：

除了调用Yalmip工具箱中的鲁棒优化模块之外，还可以使用KKT条件，以及对偶变换对该问题进行求解，也都可以调用Yalmip函数进行求解，介绍如下：

2.2 KKT条件求解

单阶段鲁棒优化问题的数学本质就是一种特殊的双层优化问题，当模型为线性或为凸优化形式时，可以通过KKT条件进行求解，KKT条件可以手动写，也可以使用Yalmip中内置的求解KKT条件函数，我在之间的博客中已有介绍(双层优化入门(2)—基于yalmip的双层优化求解(附matlab代码))，这次不再介绍，主要教大家如何手动写KKT条件并对单阶段鲁棒优化问题进行求解，相关理论如果不是很清楚的可以参考一下这篇文章(KKT条件，原来如此简单)。

为了方便求解，可以首先将上述优化问题改写成如下形式：

仔细观察这种构造的形式，因为内层优化中相当于只包含变量z，变量x为常量，因此约束条件可以只考虑含变量z的部分，其内层优化的拉格朗日函数为：

其中，λ和π是n维变量，u是1维变量。

据拉格朗日函数可以写出内层优化的KKT方程组：

将内层优化的KKT方程组添加到外层优化中，就可以将双层优化问题转为单层优化问题，如下式所示：

其中该优化问题的最后三项为互补松弛条件，是非线性约束，可以引入0-1变量，利用大M法进行等效线性化，过程如下：

其中，q,d为n维的0-1变量，v为1维的0-1变量。

简单分析一下，当qi=0时，λi = 0，zi -1 <= 0，当qi=1时，λi≥0，zi -1 = 0，满足互补松弛条件；当di=0时，πi = 0，-zi <= 0，当di=1时，πi≥0，zi = 0，满足互补松弛条件；当v=0时，u = 0， $\sum_{i=1}^nz_i-\Gamma$ ≤0，当di=1时，u≥0， $\sum_{i=1}^nz_i-\Gamma$ =0，满足互补松弛条件。因此该转换是完全等效的。至此，max-min形式的鲁棒优化问题就转变为了max形式的混合整数线性规划问题，可以快速求解。

完整的matlab求解代码如下：

%% 选股问题的鲁棒优化KKT条件求解
clc
clear
close all
yalmip('clear')

%% 1.相关的参数
n  = 150;                                               % 股票的数量 
p  = 1.15 + 0.05/150*(1:n)';                            % 期望的收益
sigma = 0.05/450*sqrt(2*n*(n + 1)*(1:n)');              % 收益的偏差
gamma = 5;                                              % 不确定预算
M = 10e2;                                               % 大M法的常数

%% 2.定义决策变量
z = sdpvar(n,1);                                        % 不确定变量z
x = sdpvar(n,1);                                        % 决策变量x
L = sdpvar(n,1);                                        % 对偶变量λ
pai = sdpvar(n,1);                                      % 对偶变量π
u = sdpvar(1);                                          % 对偶变量u
q = binvar(n,1);                                        % 决策变量q
d = binvar(n,1);                                        % 决策变量d
v = binvar(1,1);                                        % 决策变量v

%% 3.确定目标函数
objective = -(p - sigma.*z)'*x;                         % 目标函数

%% 4.确定约束条件
Constraints = [sum(x) == 1 , x >= 0];
Constraints = [Constraints , - sigma.*x + L - pai + u == 0];
Constraints = [Constraints , z >= 0 , z <= 1 , sum(z) <= gamma];
Constraints = [Constraints , L >= 0 , pai >= 0 , u >= 0];
Constraints = [Constraints , L <= M*q , z - 1 >= M*(q - 1)];
Constraints = [Constraints , pai <= M*d , -z >= M*(d - 1)];
Constraints = [Constraints , u <= M*v , sum(z) - gamma >= M*(v - 1)];
Constraints = [Constraints , boundingbox(Constraints,[],[z',x',L',pai',u])];

%% 5.求解优化问题

% gurobi求解器
% ops = sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1);
% ops.gurobi.TimeLimit = 600;                             % 运行时间限制为10min
% ops.gurobi.MIPGap = 0.024;                              % 收敛精度限制为0.024

% cplex求解器
ops = sdpsettings('verbose', 3, 'solver', 'cplex','showprogress',1);
ops.cplex.timelimit = 600;                              % 运行时间限制为10min
ops.cplex.mip.tolerances.mipgap = 0.0245;               % 收敛精度限制为0.0245

sol = optimize(Constraints,objective,ops);
if sol.problem == 0
    disp('求解成功');
else
    disp('出错了(或者超时自动结束了)');
    yalmiperror(sol.problem)
end

%% 6.优化结果
x = value(x);                                           % 决策变量x取值
plot(x)                                                 % 画出图像

运行结果如下：

可以看到和上面我们用Yalmip鲁棒优化模块的求解结果是一致的。再次总结一下使用KKT条件求解鲁棒优化问题的步骤：

1)观察鲁棒优化中内层优化和外层优化的形式，如果内层优化和外层优化目标函数一样，为了方便写出内层优化的拉格朗日函数，可以把内层优化中与内层优化变量无关的约束全部去除。具体来说，就是当内层优化的变量是不确定变量时，便只保留不确定变量相关的约束，当内层优化的变量是人工的决策变量时，便可以去除只包含不确定变量的约束。注意，如果约束条件中同时包括不确定变量和人工决策变量，则无论何种情况，该约束都需要保留。

2)引入拉格朗日乘子，写出内层优化的拉格朗日函数，并根据拉格朗日函数写出内层优化的KKT条件。这时候需要注意的是，KKT条件和优化问题取极值并不一定是完全等效的，只有当强对偶定理成立时，两者才互为充要条件。在这里教大家一个最简单的判断方式，对于常见的线性规划，二次规划与凸优化，只要原问题有解，都满足强对偶定理，则KKT条件是优化问题取极值的充要条件；而当优化问题含有0-1变量或整数变量时，原问题和对偶问题的对偶间隙不为0，此时KKT条件是优化问题取极值的必要不充分条件，如果按照相同的方式求解含0-1变量的双层优化问题，那么最终得到的结果只是一个局部最优解而不是全局最优解。因此，很多文献中的处理方式都是将优化问题分为两个阶段，将0-1变量全部包含在第一阶段中，第二阶段就是不含0-1变量的单阶段鲁棒优化问题，可以使用KKT条件方便地进行求解。

3)KKT条件中含有互补松弛条件，是一种典型的非线性约束，需要使用大M法并引入0-1变量将其转换为线性约束，转换之后可以自行验证一下当0-1变量取0或者取1时，是否和原始的互补松弛条件等效。

4)经过上述处理，可以将原始的鲁棒优化问题转为混合整数规划(如果原来是线性鲁棒优化问题，就可以转换为混合整数线性规划)，采用Yalmip+cplex，gurobi等求解器快速求解。

2.3 使用对偶变换进行求解

单阶段鲁棒优化另一种常见的求解方法是利用对偶变换，将内层优化的方向转为和外层优化相同(例如min-max问题转为min-min问题)，然后与外层优化合并，最终转为单层优化问题，并利用求解器进行求解。有关对偶转换的基本原理，我在上一篇博客中已经做了详细的讲解(双层优化入门(4)—基于对偶变换的双层优化求解)，这里不再过多介绍，再强调一下鲁棒优化可以使用对偶变换转为单层优化的条件，内层优化要为线性规划，没有非线性项及0-1变量，整数变量，否则强对偶定理不成立，对偶变换后得到的内层最优解和原内层问题的最优解不相同，就算能求出最优解，也只是一个局部最优解而不是全局最优。

同样以上面的选股优化问题为例，说明如何使用对偶变换对单阶段鲁棒优化问题进行求解，其中内层优化问题可以写作：

其中α，β是n维的对偶变量，γ是1维的对偶变量。简单分析一下对偶转换过程，原问题为min问题，所以对偶问题为max问题；原问题中有n个变量，所以对偶问题中有n个约束；原问题中有2n+1个约束条件，所以对偶问题中共有2n+1个变量，因此内层优化的对偶问题可以写作：

将内层优化的对偶问题和外层问题合并，得到

大功告成，原始的鲁棒优化问题就被转换为一个线性规划问题，用求解器就可以解了。Matlab+Yalmip求解代码如下：

%% 选股问题的鲁棒优化对偶变换求解
clc
clear
close all
yalmip('clear')

%% 1.相关的参数
n  = 150;                                               % 股票的数量 
p  = 1.15 + 0.05/150*(1:n)';                            % 期望的收益
sigma = 0.05/450*sqrt(2*n*(n + 1)*(1:n)');              % 收益的偏差
gamma = 5;                                              % 不确定预算

%% 2.定义决策变量
x = sdpvar(n,1);                                        % 决策变量x
a = sdpvar(n,1);                                        % 对偶变量α
b = sdpvar(n,1);                                        % 对偶变量β
Y = sdpvar(1);                                          % 对偶变量γ

%% 3.确定目标函数
objective = -(p'*x + sum(a) + Y*gamma);                 % 目标函数

%% 4.确定约束条件
Constraints = [sum(x) == 1 , x >= 0];
Constraints = [Constraints , a - b + Y <= -sigma.*x];
Constraints = [Constraints , a <= 0 , b <= 0 , Y <= 0];

%% 5.求解优化问题

% gurobi求解器
% ops = sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1);

% cplex求解器
ops = sdpsettings('verbose', 3, 'solver', 'cplex','showprogress',1);

sol = optimize(Constraints,objective,ops);
if sol.problem == 0
    disp('求解成功');
else
    disp('出错了(或者超时自动结束了)');
    yalmiperror(sol.problem)
end

%% 6.优化结果
x = value(x);                                           % 决策变量x取值
plot(x)                                                 % 画出图像

运行结果：

与使用鲁棒优化模块，使用KKT条件求解的结果都是一样的，说明了三种方法都是有效的。

再次总结一下使用对偶变换求解鲁棒优化问题的步骤：

1)写出鲁棒优化的内层优化问题，把内层优化中与内层优化变量无关的约束全部去除。具体来说，就是当内层优化的变量是不确定变量时，便只保留不确定变量相关的约束，当内层优化的变量是人工的决策变量时，便可以去除只包含不确定变量的约束。

2)对内层优化问题进行对偶变换，转换方式不清楚的可以参考我之前的博客(双层优化入门(4)—基于对偶变换的双层优化求解)。

3)将内层优化问题和外层优化问题合并，得到一个单层优化问题，可采用Gurobi等求解器进行求解。

3.电力系统鲁棒经济调度Matlab+Yalmip编程实例

参考文献[1]中将单阶段鲁棒优化问题描述为人工决策与大自然的博弈问题，这么做在说法上具有一定的创新性，但问题的本质还是一样，万变不离其宗，只要掌握了求解方法，不管鲁棒优化的问题背景怎么样，只要可以转为相同的形式，就可以采用相同的方法进行求解。

文献[1]中所述的含风电和电动汽车的电力系统鲁棒经济调度数学模型为：

我们还是和鲁棒选股优化问题一样，使用三种不同的方式求解这个鲁棒优化问题，分别如下。

3.1 鲁棒优化模块求解

这是一个实际的优化问题，看起来相当复杂，但是也不用慌，我们理清思路，分分钟拿下。

1）优化问题中的决策变量有哪些？

2）优化问题中的不确定变量是什么？

不确定变量是风电场k在t时刻的可用风功率 $p_{wk,t}^{w}$ ，是一个Nw×T维的变量。此外，对于比较复杂的问题，我在Matlab编程之前通常都会将问题中的参数和变量分别用表格罗列出来，并写出各个变量和参数在Matlab的定义。这个鲁棒优化问题的参数和变量分别如表1和表2所示：

表1 决策变量

表2 相关参数

另外，为保证系统潮流收敛，变量pl,t需要添加到决策变量中，同时还要列写潮流平衡方程。因为文献中节点编号和MATPOWER工具箱中的编号有些不一样，为了编程方便，我统一改为MATPOWER中的编号形式，并规定了潮流方向，如下图所示：

由于该系统只有9个节点，因此可以把潮流方程直接列出来，具体如下：

3）确定优化问题中不确定集合的形式，并考虑是否可以直接使用Yalmip中的函数进行求解。

优化问题中不确定集合为简单的盒式不确定集，可以直接使用Yalmip的鲁棒优化功能进行求解。

4）确定优化问题的目标函数。

5）确定优化问题中的约束条件。

经过上面5个步骤的处理，便顺利完成了数学建模，并可以使用Yalmip中的鲁棒优化模块求解这个问题，代码在压缩包中的Problem2文件夹中，运行Problem2_robust.m文件即可得到结果。

运行结果如下：

由结果可知，大自然的风功率在所有时刻都取得下限值，也就是鲁棒优化中最恶劣的场景，此时大电网调度的风功率和大自然的风功率完全相等，说明这个方法在最恶劣的场景下也能得到最优的调度方案(由于Yalmip工具箱中的鲁棒优化模块求解时并不会返回不确定变量的取值，所以没有画出大自然的风功率，但是由弃风和备用费用之和为0可知大自然的风功率和大电网调度的风功率时相等的)。

3.2 KKT条件求解

我们首先写出上面鲁棒优化问题的内层优化问题，由于内层优化是和不确定变量 $p_{wk,t}^{W}$ 有关，我们可以只列出和该变量相关的目标函数及约束，具体如下：

这个优化问题仅含有变量 $p_{wk,t}^{W}$ ，相对比较简单，但是目标函数中包含非线性函数max，首先需要进行等效线性化，我们引入两个中间变量 $p_{wk,t}^{W+}$ 和 $p_{wk,t}^{W-}$ ，分别表示变量 $p_{wk,t}^{W}$ 中大于和小于 $p_{wk,t}^{G}$ 的分量，即 $p_{wk,t}^{W+}=\max(0,p_{wk,t}^{W}-p_{wk,t}^{G})$ 以及 $p_{wk,t}^{W-}=\max(0,p_{wk,t}^{G}-p_{wk,t}^{W})$ 。

注意，和选股鲁棒优化问题不同的是，这个内层优化是一个max问题，因此需要把所有的约束都写成≥0的形式，再求KKT条件。转换后没有非线性项，可以使用KKT条件进行等效转换。我们首先写出拉格朗日函数：

其中， $\alpha,\beta,\gamma,\lambda,\omega,\pi$ 是拉格朗日乘子，都是Nw×T维的变量。进一步写出内层优化问题的KKT条件：

另外，KKT条件中引入了互补松弛条件，为非线性约束，需要通过引入0-1变量进行转换，转换思路和2.2节一样，这里不再赘述，只给出最终转换结果：

其中， $q,d,u,\nu,\eta,\theta$ 均为Nw×T维的0-1变量，这样处理，就把内层优化问题处理成了无目标的KKT条件，将其添加到上层优化问题的目标函数中，就可以将原来的鲁棒优化问题转换为混合整数二次规划问题：

使用Yalmip+Gurobi求解上述混合整数线性规划问题的代码在压缩包中的Problem2文件夹中，运行Problem2_kkt.m文件即可得到结果，运行结果如下：

好家伙，得到的结果居然和Yalmip鲁棒优化模块直接求解的结果不一样，而且此时大自然决定的风功率都取得各个区间的上限值，这不是最恶劣的场景，而是最乐观的场景好吧.......我仔细检查模型很多遍，也并没有发现问题出在哪里，也许可能是某个线性转换的过程存在误差吧。寻找问题的过程也是对巩固理论知识的一种方式，大家可以帮忙看看，如果发现问题的话私信我一下。

还有一件诡异的事情，调试代码的时候我有一次把中间变量 $p_{wk,t}^{W+}$ 和 $p_{wk,t}^{W-}$ 的上下限约束删除了，然后运行得到了一个无界解，居然和文献[1]中不确定变量的结果有点相似，也就是自然界的风功率在上下限之间反复横跳(不太理解为什么求出的最恶劣场景会在上下限之间反复横跳，按我常规的脑回路来考虑，都取下限值时应该就是最恶劣的场景)。

当然，这很明显不是最优解，贴出来仅供娱乐。错误的代码我也放在了Problem2_kkt1.m文件中，运行即可得到结果。

3.3 使用对偶变换进行求解

上面用KKT条件求解时已经写出了内层优化的线性化表达，直接拿过来，并稍做变换：

其中α—η都是对偶变量，对偶转换过程如下，

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

2)原问题中有20个变量(把Nw=1，T=10看成已知数)，因此对偶问题中有20个约束条件。

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

4)原问题中目标函数的系数分别为10个Cu,wk和10个Co,wk,因此对偶问题的约束条件中≥号右边的常数也为10个Cu,wk和10个Co,wk。

5)原问题中有60个约束条件，因此对偶问题有60个决策变量，即α—η都是10维的决策变量。

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

7)对偶问题的约束条件中各个对偶变量的系数和原问题约束条件中变量的系数互为转置关系。

综上所述，内层优化的对偶问题可以写作：

将内层优化的对偶问题和外层问题合并，得到

经过对偶变换后，min-max类型的鲁棒优化问题被转为单层的min问题，但不幸的是，目标函数也从线性转为非线性，变成了一个二次函数，而这个问题约束条件均为线性约束，这是1个二次规划问题。幸运的是，二次规划是最简单的一种非线性规划，求解方式有很多，比如可以使用KKT条件进行求解(感兴趣的朋友可以尝试写一下)，而且外层优化的目标函数本身包含一个二次函数，因此并没有给模型的求解增加任何困难。当然最便捷的方式还是调用Cplex、Gurobi等求解器，在这里我直接使用求解器进行求解。

使用Yalmip+Gurobi求解上述二次规划问题的代码在压缩包中的Problem2文件夹中，运行Problem2_dual.m文件即可得到结果。具体运行结果如下：

通过对偶变换求解鲁棒优化问题的结果和KKT条件求解结果一样，而且此时大自然决定的风功率都是取得各个区间的上限值，不像是最恶劣场景，而像是最乐观的场景，和我们的常理推断，以及使用鲁棒优化模块直接求解的结果都不一样，到底哪个是最优解。。。。说实话我也搞不清楚了，大家自行判断吧。

4.总结

通过两个单阶段鲁棒优化的实例，系统介绍了如何通过Matlab+Yalmip工具箱求解单阶段鲁棒优化的方法，总共包括三种方法，一是使用Yalmip中的uncertain函数定义不确定变量，并直接使用鲁棒优化模块进行求解；二是利用将内层优化取最优解的KKT条件到外层优化中，转为单层优化进行求解，三是利用对偶变换，将内外层优化的目标函数的方向调整为一致，合并形成单层优化。其中第1个简单的鲁棒选股优化问题(看起来简单，但其实变量规模更大)，三种方法优化结果一致，而第2个电力系统鲁棒经济调度问题，KKT条件和对偶变换的求解结果一致，但和直接调用鲁棒优化模块的结果不同，具体原因还在探索之中。三种方法都有各自的适用条件：

1.想要用uncertain函数定义不确定变量，并通过鲁棒优化模块直接求解，不确定集合需满足一定要求。可以是盒式不确定集、多面体不确定集、椭球不确定集以及1，2，∞范数约束的不确定集合。另外，似乎对于数据驱动的分布式鲁棒优化，也可以使用Yalmip进行求解(见官方文档Sample-based robust optimization - YALMIP)。但我现在还没研究，有空我学会了再来和大家分享经验。

2.对于KKT条件和对偶转换的方法，两者都是对模型有一定的要求，需要内层优化问题满足强对偶定理(只要包含0-1变量或整数变量，一定不满足，即不能通过KKT条件和对偶变换求最优解)。要想准确判断，需要一定的数学基础，对于实际应用来说，只要知道当内层优化问题的约束条件都为线性，目标函数为线性函数或者二次函数，且决策变量不包含0-1变量和整数型变量时，一般都可以通过KKT条件和对偶转换进行求解。