鲁棒优化入门（二）——基于matlab+yalmip求解鲁棒优化问题

       上一篇博客简单介绍了可以用来求解鲁棒优化的两个工具箱：

鲁棒优化入门（一）——工具箱Xprog和RSOME的安装与使用

        其实大家可能没有想过，matlab+yalmip工具箱也可以处理一些简单的鲁棒优化问题，上官方文档：Robust optimization - YALMIP

        这里就和大家一起学习一下使用yalmip+cplex求解鲁棒优化问题的方法。

一、yalmip求解鲁棒优化

        鲁棒优化问题可以表示为如下的一般形式：

​

式中，x表示决策变量，w表示不确定变量。这个表达式的含义就是在最恶劣的情况下（w的取值使目标函数最大），求出满足约束条件，并且能使目标函数最小的决策变量x。

        但yalmip工具箱求解鲁棒优化问题时还是有一定局限性的，并不能处理任意形式的不确定集下的鲁棒优化问题，一般来说，当鲁棒优化问题的不确定集合为箱型不确定集、椭球不确定集以及多面体不确定集时，都可以用yalmip工具箱求解（具体细节可参考官方文档）。

二、示例代码

1.实例1：线性规划问题

        考虑一个简单线性规划问题：

​

        这个问题中不确定变量w仅存在于约束条件中，目标函数中不包含不确定变量。显然，当w=0.5时是最恶劣的情况，此时f(x)的最大值为0.5。采用yalmip编程验证一下：

sdpvar x w                              % 定义变量
C = [x+w <= 1];                         % 约束条件
W = [-0.5 <= w <= 0.5, uncertain(w)];   % 不确定集
objective = -x;                         % 目标函数
sol = optimize(C + W,objective);        % 求解模型
obj=-value(objective);                  % 目标函数取值
x=value(x);                             % 决策变量x取值

        对于不确定变量w，需要用uncertain()函数将其规定为不确定变量。另外，还可以先将鲁棒优化存为yalmip模型，然后再进行求解：

[Frobust,robust_objective] = robustify(C + W,objective);    % 导出鲁棒优化模型
sol = optimize(Frobust,robust_objective);                   % 求解鲁棒优化模型

运行结果：

         显然，和我们一眼看出的结果是一样的。

        对于决策变量为整数或含有逻辑约束的鲁棒优化问题，yalmip同样可以求解，例如：

​

         matlab代码为：

intvar x
sdpvar w                                                    % 定义变量
C = [x+w <= 2];                                             % 约束条件
W = [-0.5 <= w <= 0.5, uncertain(w)];                       % 不确定集
objective = -x;                                             % 目标函数
sol = optimize(C + W,objective);                            % 求解模型
obj=-value(objective);                                      % 目标函数取值
x=value(x);                                                 % 决策变量x取值

        运行结果：

2.实例2：含椭球不确定集的鲁棒优化问题

        考虑一个含有椭球不确定集的鲁棒优化问题：

​

         这个问题中不确定变量w仅存在于约束条件中，目标函数中不包含不确定变量。假设n=2，也不难看出，当w1=w2=0.5时，属于最恶劣的场景，此时f(x)最大值为0。编程验证一下：

sdpvar x w(2,1)                                             % 定义变量
C = [x+sum(w) <= 1];                                        % 约束条件
W = [norm(w) <= 1/sqrt(2), uncertain(w)];                   % 不确定集
objective = -x;                                             % 目标函数
sol = optimize(C + W,objective);                            % 求解模型
obj=-value(objective);                                      % 目标函数取值
x=value(x);                                                 % 决策变量x取值

        运行结果：

        小伙伴们可能有点奇怪，说好的最优值是零呢？怎么是一个负数？其实，我们应该知道matlab的计算精度是有限的，eps表示MATLAB默认的最小浮点数精度(默认是eps(1))：

​

         所以这个结果应该就是计算误差。

3.实例3：含不确定性的平方和(SOS,Sum of squares)规划问题

        考虑如下的SOS规划问题，其中a为整数变量，取值范围[3,5]

​

         这玩意就没法直接看出最优解了，直接上代码：

intvar a
sdpvar x y t u                  % 定义决策变量
p = a*x^4+y^4+u*x*y+1;          % 多项式表示
F = [uncertain(u), -1<=u<=1];   % 不确定集
F = [F, a>=3, a<=5];            % 约束条件a∈[3,5]
F = [F, sos(p-t)];              % 多项式约束
solvesos(F, -t)                 % 求解模型
A=value(a);                     % 决策变量a的取值 
T=value(t);                     % 目标函数t的取值

运行结果：

        可以看出，当a=5时，目标函数t取得最大值，为0.9437。

        上面三个例子都是yalmip官方文档中的示例，下面来看一个实际问题：

4.实例4：股票投资问题

        假设一共有150种股票可供选择，第i个股票的不确定收益用表示，其取值满足约束条件：，其中，表示股票i的期望收益，表示股票i的偏差，不确定集可以用1范数和无穷范数表示为：

         该投资组合问题的鲁棒优化模型可以表示为：

         在本例中，假设Γ=5，参数pi和σi满足如下公式：

         求解该问题的matlab代码如下：

n  = 150;                                               % 股票的数量 
p  = 1.15+ 0.05/150*(1:n)';                             % 期望的收益
sigma = 0.05/450*sqrt(2*n*(n+1)*(1:n)');                % 收益的偏差
gamma=5;                                                % 不确定预算
z=sdpvar(n,1);                                          % 不确定变量z
x=sdpvar(n,1);                                          % 决策变量x
C=[sum(x)==1,x>0];                                      % x的约束条件
Z=[norm(z,Inf)<=1,norm(z,1)<=gamma,uncertain(z)];       % 不确定集
objective = -(p + sigma.*z)'*x;                         % 目标函数
sol = optimize(C+Z,objective);                          % 求解模型
x=value(x);                                             % 决策变量x取值
plot(x)                                                 % 画出图像

        运行结果：

        和RSOME工具箱的求解结果一致。