商场网站模板,免费做金融网站,建一个个人网站要多少钱,好听的公司名字之前的博客#xff1a;鲁棒优化入门#xff08;二#xff09;——基于matlabyalmip求解鲁棒优化问题 去年发布了使用Yalmip工具箱求解鲁棒优化问题的博客之后#xff0c;陆陆续续有朋友问我相关的问题#xff0c;有人形容从学习这篇博客到求解论文中的鲁棒优化问题#x… 之前的博客鲁棒优化入门二——基于matlabyalmip求解鲁棒优化问题 去年发布了使用Yalmip工具箱求解鲁棒优化问题的博客之后陆陆续续有朋友问我相关的问题有人形容从学习这篇博客到求解论文中的鲁棒优化问题就好像刚学会求导公式就要去做高考压轴题根本无从下手。为了解决这个问题这篇博客将手把手地教会大家如何使用Matlab yalmipcplex(当然其他的求解器比如gurobi也是可以的)求解论文中的鲁棒优化问题。为了具有拓展性本文将选取两篇不同的单阶段鲁棒优化问题一个是经典的选股优化问题另一个是参考文献[1]中所提的电力系统鲁棒经济调度问题分别说明如何对鲁棒优化的文献进行分析并借助Yalmip工具箱实现文献中鲁棒优化问题的编程求解。
1.单阶段鲁棒优化基本形式 如参考文献[1]中所述标准的单阶段鲁棒优化的形式为 其中x表示系统的状态(也可以认为是参数)u表示人为控制变量(也就是决策变量)w表示大自然的决策变量(也就是不确定变量)。万变不离其宗几乎所有的单阶段鲁棒优化问题都可以写成这种形式并采用相同的套路进行求解。因此对于任意文献中的单阶段鲁棒优化问题可以转为和上式统一的形式便于求解具体转换思路如下 1确定优化问题中的决策变量有哪些将其与变量u对应。 2确定优化问题中的不确定变量是什么将其与变量w对应。 3确定优化问题中不确定集合的形式并考虑是否可以直接使用Yalmip中的鲁棒优化模块进行求解。 4确定优化问题的目标函数将其与J(u,w)对应。 5确定优化问题中的约束条件是否包含非线性约束需要线性化。
2.选股优化问题MatlabYalmip编程实例 以一个经典的选股优化问题为例说明如何使用Matlab实现编程求解单阶段鲁棒优化问题。 一共有n150支股票其中第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变量。 简单分析一下当qi0时λi 0zi -1 0当qi1时λi≥0zi -1 0满足互补松弛条件当di0时πi 0-zi 0当di1时πi≥0zi 0满足互补松弛条件当v0时u 0≤0当di1时u≥00满足互补松弛条件。因此该转换是完全等效的。至此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.0245sol optimize(Constraints,objective,ops);
if sol.problem 0disp(求解成功);
elsedisp(出错了(或者超时自动结束了));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)经过上述处理可以将原始的鲁棒优化问题转为混合整数规划(如果原来是线性鲁棒优化问题就可以转换为混合整数线性规划)采用Yalmipcplexgurobi等求解器快速求解。
2.3 使用对偶变换进行求解 单阶段鲁棒优化另一种常见的求解方法是利用对偶变换将内层优化的方向转为和外层优化相同(例如min-max问题转为min-min问题)然后与外层优化合并最终转为单层优化问题并利用求解器进行求解。有关对偶转换的基本原理我在上一篇博客中已经做了详细的讲解(双层优化入门(4)—基于对偶变换的双层优化求解)这里不再过多介绍再强调一下鲁棒优化可以使用对偶变换转为单层优化的条件内层优化要为线性规划没有非线性项及0-1变量整数变量否则强对偶定理不成立对偶变换后得到的内层最优解和原内层问题的最优解不相同就算能求出最优解也只是一个局部最优解而不是全局最优。 同样以上面的选股优化问题为例说明如何使用对偶变换对单阶段鲁棒优化问题进行求解其中内层优化问题可以写作 其中αβ是n维的对偶变量γ是1维的对偶变量。简单分析一下对偶转换过程原问题为min问题所以对偶问题为max问题原问题中有n个变量所以对偶问题中有n个约束原问题中有2n1个约束条件所以对偶问题中共有2n1个变量因此内层优化的对偶问题可以写作 将内层优化的对偶问题和外层问题合并得到 大功告成原始的鲁棒优化问题就被转换为一个线性规划问题用求解器就可以解了。MatlabYalmip求解代码如下
%% 选股问题的鲁棒优化对偶变换求解
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 0disp(求解成功);
elsedisp(出错了(或者超时自动结束了));yalmiperror(sol.problem)
end%% 6.优化结果
x value(x); % 决策变量x取值
plot(x) % 画出图像 运行结果 与使用鲁棒优化模块使用KKT条件求解的结果都是一样的说明了三种方法都是有效的。 再次总结一下使用对偶变换求解鲁棒优化问题的步骤 1)写出鲁棒优化的内层优化问题把内层优化中与内层优化变量无关的约束全部去除。具体来说就是当内层优化的变量是不确定变量时便只保留不确定变量相关的约束当内层优化的变量是人工的决策变量时便可以去除只包含不确定变量的约束。 2)对内层优化问题进行对偶变换转换方式不清楚的可以参考我之前的博客(双层优化入门(4)—基于对偶变换的双层优化求解)。 3)将内层优化问题和外层优化问题合并得到一个单层优化问题可采用Gurobi等求解器进行求解。
3.电力系统鲁棒经济调度MatlabYalmip编程实例 参考文献[1]中将单阶段鲁棒优化问题描述为人工决策与大自然的博弈问题这么做在说法上具有一定的创新性但问题的本质还是一样万变不离其宗只要掌握了求解方法不管鲁棒优化的问题背景怎么样只要可以转为相同的形式就可以采用相同的方法进行求解。 文献[1]中所述的含风电和电动汽车的电力系统鲁棒经济调度数学模型为 我们还是和鲁棒选股优化问题一样使用三种不同的方式求解这个鲁棒优化问题分别如下。
3.1 鲁棒优化模块求解 这是一个实际的优化问题看起来相当复杂但是也不用慌我们理清思路分分钟拿下。 1优化问题中的决策变量有哪些 2优化问题中的不确定变量是什么 不确定变量是风电场k在t时刻的可用风功率是一个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条件求解 我们首先写出上面鲁棒优化问题的内层优化问题由于内层优化是和不确定变量有关我们可以只列出和该变量相关的目标函数及约束具体如下 这个优化问题仅含有变量 相对比较简单但是目标函数中包含非线性函数max首先需要进行等效线性化我们引入两个中间变量和分别表示变量中大于和小于的分量即以及。 注意和选股鲁棒优化问题不同的是这个内层优化是一个max问题因此需要把所有的约束都写成≥0的形式再求KKT条件。转换后没有非线性项可以使用KKT条件进行等效转换。我们首先写出拉格朗日函数 其中是拉格朗日乘子都是Nw×T维的变量。进一步写出内层优化问题的KKT条件 另外KKT条件中引入了互补松弛条件为非线性约束需要通过引入0-1变量进行转换转换思路和2.2节一样这里不再赘述只给出最终转换结果 其中 均为Nw×T维的0-1变量这样处理就把内层优化问题处理成了无目标的KKT条件将其添加到上层优化问题的目标函数中就可以将原来的鲁棒优化问题转换为混合整数二次规划问题 使用YalmipGurobi求解上述混合整数线性规划问题的代码在压缩包中的Problem2文件夹中运行Problem2_kkt.m文件即可得到结果运行结果如下 好家伙得到的结果居然和Yalmip鲁棒优化模块直接求解的结果不一样而且此时大自然决定的风功率都取得各个区间的上限值这不是最恶劣的场景而是最乐观的场景好吧.......我仔细检查模型很多遍也并没有发现问题出在哪里也许可能是某个线性转换的过程存在误差吧。寻找问题的过程也是对巩固理论知识的一种方式大家可以帮忙看看如果发现问题的话私信我一下。 还有一件诡异的事情调试代码的时候我有一次把中间变量和的上下限约束删除了然后运行得到了一个无界解居然和文献[1]中不确定变量的结果有点相似也就是自然界的风功率在上下限之间反复横跳(不太理解为什么求出的最恶劣场景会在上下限之间反复横跳按我常规的脑回路来考虑都取下限值时应该就是最恶劣的场景)。 当然这很明显不是最优解贴出来仅供娱乐。错误的代码我也放在了Problem2_kkt1.m文件中运行即可得到结果。
3.3 使用对偶变换进行求解 上面用KKT条件求解时已经写出了内层优化的线性化表达直接拿过来并稍做变换 其中α—η都是对偶变量对偶转换过程如下 1)原问题为max问题对偶问题为min问题。 2)原问题中有20个变量(把Nw1T10看成已知数)因此对偶问题中有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等求解器在这里我直接使用求解器进行求解。 使用YalmipGurobi求解上述二次规划问题的代码在压缩包中的Problem2文件夹中运行Problem2_dual.m文件即可得到结果。具体运行结果如下 通过对偶变换求解鲁棒优化问题的结果和KKT条件求解结果一样而且此时大自然决定的风功率都是取得各个区间的上限值不像是最恶劣场景而像是最乐观的场景和我们的常理推断以及使用鲁棒优化模块直接求解的结果都不一样到底哪个是最优解。。。。说实话我也搞不清楚了大家自行判断吧。
4.总结 通过两个单阶段鲁棒优化的实例系统介绍了如何通过MatlabYalmip工具箱求解单阶段鲁棒优化的方法总共包括三种方法一是使用Yalmip中的uncertain函数定义不确定变量并直接使用鲁棒优化模块进行求解二是利用将内层优化取最优解的KKT条件到外层优化中转为单层优化进行求解三是利用对偶变换将内外层优化的目标函数的方向调整为一致合并形成单层优化。其中第1个简单的鲁棒选股优化问题(看起来简单但其实变量规模更大)三种方法优化结果一致而第2个电力系统鲁棒经济调度问题KKT条件和对偶变换的求解结果一致但和直接调用鲁棒优化模块的结果不同具体原因还在探索之中。三种方法都有各自的适用条件 1.想要用uncertain函数定义不确定变量并通过鲁棒优化模块直接求解不确定集合需满足一定要求。可以是盒式不确定集、多面体不确定集、椭球不确定集以及12∞范数约束的不确定集合。另外似乎对于数据驱动的分布式鲁棒优化也可以使用Yalmip进行求解(见官方文档Sample-based robust optimization - YALMIP)。但我现在还没研究有空我学会了再来和大家分享经验。 2.对于KKT条件和对偶转换的方法两者都是对模型有一定的要求需要内层优化问题满足强对偶定理(只要包含0-1变量或整数变量一定不满足即不能通过KKT条件和对偶变换求最优解)。要想准确判断需要一定的数学基础对于实际应用来说只要知道当内层优化问题的约束条件都为线性目标函数为线性函数或者二次函数且决策变量不包含0-1变量和整数型变量时一般都可以通过KKT条件和对偶转换进行求解。
5.完整代码 代码使用matlabYalmip求解第二个鲁棒优化问题中还用到了MATPOWER工具箱下载链接如下
鲁棒优化入门(5)-MatlabYalmip求解鲁棒优化编程实战
参考文献
[1]梅生伟,郭文涛,王莹莹等.一类电力系统鲁棒优化问题的博弈模型及应用实例[J].中国电机工程学报,2013,33(19):47-5620.
[2]KKT条件原来如此简单
[3]鲁棒优化入门(2)—基于MatlabYalmip求解鲁棒优化
[4]双层优化入门(2)—基于yalmip的双层优化求解(附matlab代码)
[5]双层优化入门(4)—基于对偶变换的双层优化求解
[6]Sample-based robust optimization - YALMIP
