肇庆住房和城乡建设局网站,龙华网站建设的软件,免费的建网站软件,广告推广赚钱在哪接0.引言 上一篇博客介绍了使用Yalmip工具箱求解单阶段鲁棒优化的方法。这篇文章将和大家一起继续研究如何使用Yalmip工具箱求解两阶段鲁棒优化(默认看到这篇博客时已经有一定的基础了#xff0c;如果没有可以看看我专栏里的其他文章)。关于两阶段鲁棒优化与列与约束生成算法的原…0.引言 上一篇博客介绍了使用Yalmip工具箱求解单阶段鲁棒优化的方法。这篇文章将和大家一起继续研究如何使用Yalmip工具箱求解两阶段鲁棒优化(默认看到这篇博客时已经有一定的基础了如果没有可以看看我专栏里的其他文章)。关于两阶段鲁棒优化与列与约束生成算法的原理之前的博客已经详细地介绍过了这里就不再过多介绍主要是结合实例来讲解编程思路。这篇博客用到了两个算例1个是两阶段鲁棒优化问题和列与约束生成算法的开山鼻祖[1]另一个是电气专业中两阶段鲁棒优化问题最热门的文章之一[2]相信大家在网上见到过无数号称完美复现的代码但实际上大部分都是有问题的(包括我自己早期写的代码也是被网上的代码带歪了后面理解慢慢深入才发现问题所在)。 求解两阶段鲁棒优化问题一共有两个难点一是求解max-min或者min-max形式的子问题其实就是求解一个单阶段鲁棒优化上一篇博客我已经非常详细地介绍了求解方式借助Yalmip工具箱共有三种不同的方式可以解决。二是主问题和子问题的迭代求解也就是列与约束生成算法(CCG)的实现。很多代码在复现CCG算法时并没有向主问题同时添加列(变量)和约束这也是代码中最常见的问题。针对这两个难点我将用两个不同的算例详细地进行讲解。 此外文献[1]和[2]中都是采用了先将约束条件写成紧凑的矩阵形式然后再对子问题进行处理的方式很多朋友和我反映这部分太难处理了实际问题的约束建模过程中经常包括循环语句想要转成矩阵形式确实很不容易。这篇文章中我将分别采用两种不同的方式求解鲁棒优化。一是采用原始的约束条件省去将约束条件转为矩阵形式的步骤这种方式数学公式可能会更繁琐但比起矩阵形式的转换理解起来会更容易一些。二是采用矩阵形式进行编程在博客中我教大家一种非常简单就能将约束写为矩阵形式的方法文中只是介绍了如何使用之后也会单独写博客对此详细展开。 总之这篇博客干货满满可以认真通读一遍跟着博客中的思路亲自动手使用Matlab Yalmip实现两阶段鲁棒优化的编程(博客中提到的所有例子我都提供了相应的代码)。相信大家理解后面对任何类型的两阶段鲁棒优化问题都能迅速使用类似的方法进行解决。 博客中主要包含8大内容 ①.拿到一个复杂的两阶段鲁棒优化问题的分析步骤和方法。 ②.采用Yalmip工具箱中的uncertain函数和鲁棒优化模块求解两阶段鲁棒优化的子问题。 ③.Yalmip工具箱中的鲁棒优化模块和常规的求解思路有什么异同。 ④.使用KKT条件求解两阶段鲁棒优化的子问题并使用CCG算法进行迭代求解。 ⑤.使用对偶变换求解两阶段鲁棒优化的子问题并使用CCG算法进行迭代求解。 ⑥.采用Yalmip工具箱的内置函数将线性约束写成紧凑矩阵形式的方法。 ⑦.矩阵形式的两阶段鲁棒优化问题如何快速写出子问题内层优化的KKT条件并使用CCG算法进行迭代求解。 ⑧.矩阵形式的两阶段鲁棒优化问题如何快速写出子问题内层优化的对偶问题并使用CCG算法进行迭代求解。 由于博客篇幅较长将分上下两篇发布其中上篇使用的是文献[1]中的算例包含上述①-⑤的内容。下篇使用文献[2]中的算例包含上述①、④、⑥-⑧的内容。 这篇博客是上篇的内容。
1.两阶段鲁棒优化基本形式 如文献[1]中所述标准的两阶段鲁棒优化问题的形式为 其中y为第一阶段决策变量u为不确定变量x为第二阶段决策变量。和分析单阶段鲁棒优化问题的五个特征一样拿到一个复杂的两阶段鲁棒优化问题先不用慌按照下面的步骤进行分析即可 1确定第一阶段决策变量有哪些将其与变量y对应。 2确定第二阶段决策变量有哪些将其与变量x对应。 3确定不确定变量有哪些将其与变量u对应。 4确定优化问题中不确定集合的形式并考虑是否可以直接使用Yalmip中的鲁棒优化模块进行求解。 5确定目标函数是否有仅包含第一阶段决策变量的项如果有的话可以单独拿出来。 6确定子问题的目标函数将其与鲁棒优化的标准形式相对应。 7确定约束条件考虑是否包含非线性约束是否需要线性化。 8求解max-min或者min-max类型的子问题。 9使用迭代方式将子问题产生的变量和约束不断添加到主问题中最终得到最优解。 下面分别以文献[1]和[2]中的优化问题进行讲解说明
2.两阶段鲁棒运输问题编程实战 文献[1]中算例分析部分运输问题的两阶段鲁棒优化模型如下 我之前写过一篇博客解析这篇论文但是使用的是紧凑的矩阵形式。由于该问题比较简单紧凑形式和一般约束形式差别不大这次博客将使用一般形式的约束进行求解方便大家体会两阶段鲁棒优化原理。
2.1 鲁棒优化模块求解 按上面的思路逐步进行分析 1确定第一阶段决策变量有哪些。 第一阶段决策变量为y和z其中yi是m维的0-1变量取1时表示在i地建造仓库z是m维的连续变量表示仓库的容量(m是待选址仓库的数目)。 2确定第二阶段决策变量有哪些。 第二阶段决策变量为xij是一个m×n的连续变量表示从i仓库运往j用户的商品数目(n是用户数)。 3确定不确定变量有哪些。 不确定变量为dj表示j用户的不确定商品需求。 4确定优化问题中不确定集合的形式并考虑是否可以直接使用Yalmip中的鲁棒优化模块进行求解。 该优化问题的不确定集中为多面体形式可以直接使用Yalmip中的鲁棒优化模块进行求解。 5确定目标函数是否有仅包含第一阶段决策变量的项如果有的话可以单独拿出来。 目标函数中有一部分仅包含第一阶段的决策变量第二阶段子问题的目标函数可以不考虑这部分。 6确定子问题的目标函数将其与两阶段鲁棒优化的标准形式相对应。 其中子问题的目标函数为: 7确定约束条件考虑是否包含非线性约束是否需要线性化。 子问题中的约束条件均为线性且决策变量中不包含0-1变量满足强对偶定理KKT条件和对偶变换都是适用的。可以直接通过uncertain函数直接使用Yalmip鲁棒优化模块。 分析完成后下面可以开始尝试求解子问题相关参数如下 需要注意的是为了子问题的模型可行需要保证三个仓库容量的总和大于用户最高的需求(也就是当g0g1g21.8时的用户需求)也就是需要添加约束条件 根据上面的公式我们可以写出各个参数矩阵以及变量的表达式 用matlab代码表示如下
%% 参数矩阵
f [400; 414; 326];
a [18; 25; 20];
k 800;
C [22, 33, 24;33, 23, 30;20, 25, 27];
d_ [206; 274; 220];
d_wave 40;
gamma [1.8,1.2];
P [1 1;1 1;1 0];%% 决策变量
y binvar(3,1);
z sdpvar(3,1);
x sdpvar(3,3,’full’);
d sdpvar(3,1);
g sdpvar(3,1);%% 目标函数
objective f*y a*z sum(sum(C.*x));%% 约束条件
Constraints [];
Constraints [Constraints , z 0 , x 0 , g 0 , g 1];
Constraints [Constraints , z k*y , sum(z) sum(d_) gamma(1)*d_wave];
Constraints [Constraints , sum(x) z];
Constraints [Constraints ,sum(x,2) d];
Constraints [Constraints ,d d_ g*d_wave];
Constraints [Constraints ,g*P gamma];%% 设置求解器
ops sdpsettings(verbose, 3, solver, gurobi);
sol optimize(Constraints,objective,ops); 可以先尝试求解一下确定性优化问题和后面的两阶段鲁棒优化进行对比 8求解max-min或者min-max类型的子问题。 为了便于调试我们首先把子问题给解决了再通过迭代求解两阶段鲁棒优化问题。其中在子问题中第一阶段的变量y和z实际都是已知量此时子问题可以转为 注意我在模型中通过等式约束消去了决策变量d并使用决策变量g来表示不确定性减少变量的数目加快求解效率。 此外由于子问题中可以将主问题的决策变量视为常数因此只含有变量yz的约束都可以省略如果变量x的约束中带有变量yz那么视为常数即可。为了子问题调试方便我们先把确定性优化结果中z的取值为[772;0;0]并带入求解子问题求解成功了再和主问题进行交互迭代Matlab代码如下
%% 两阶段鲁棒优化的子问题—采用鲁棒优化模块求解
clc
clear
close all
warning off%% 参数设置
f [400; 414; 326]; % 仓库建设费用
a [18; 25; 20]; % 单位容量存储费用
C [22, 33, 24; % 从仓库i到用户j的单位运输费用33, 23, 30;20, 25, 27];
d_ [206; 274; 220]; % 基准用户需求
d_wave 40; % 需求波动
gamma [1.8,1.2]; % 不确定预算
P [1 1;1 1;1 0]; % 不确定集合的系数矩阵%% 决策变量
z [772;0;0]; % 仓库的容量
x sdpvar(3,3,full); % 从仓库i到用户j运输的商品数
g sdpvar(3,1); % 用户需求的波动幅度%% Objective
objective sum(sum(C.*x));%% 约束条件
G [uncertain(g) , g*P gamma , g 0 , g 1 ];
Constraints [];
Constraints [Constraints , sum(x,2) z ,x 0 ];
Constraints [Constraints , sum(x) d_ g*d_wave];%% 求解优化问题
opssdpsettings(verbose, 3, solver, cplex );
soloptimize(Constraints G , objective ,ops);%% 判断求解是否成功
if sol.problem 0disp(求解成功!!!);
elsedisp([求解失败原因为,sol.info]);
end%% 优化结果
if sol.problem 0objective value(objective)y [1;0;1]x value(x)z value(z)F f*y a*z sum(sum(C.*x))
end 运行结果如下 结果显示优化问题不可行无法使用Yalmip的鲁棒优化模块进行求解。但是按照官方文档说法是可以使用鲁棒优化模块求解的针对这个问题我去询问了YALMIP工具箱的作者Johan Löfberg教授他的回复是这样的 按照老师的说法Yalmip中的鲁棒优化模块考虑的并不是求出一种最恶劣场景下的决策方案而是求出所有可能最恶劣的场景下都能满足约束条件的一种决策方案。 老师给了一个例子
x w,y 1-w,xy0.5 其中w是不确定变量对于任意一种相对恶劣的场景这个问题是具有可行解的。例如w0w1但是鲁棒优化模块求解的方案要求对所有可能最恶劣的场景决策变量都要满足约束条件因此w0时的约束条件(x0,y1)与w1时的约束条件(x1,y0)都要满足也就是x0,y0和另一个约束条件xy0.5互相冲突因此这个简单的鲁棒优化实例是没有可行解的。编程验证一下
%% Johan Löfberg教授提供的例子
clc
clearsdpvar x y wobjective x y;Constraints [uncertain(w) , w 0 , w 1];
Constraints [Constraints , x w , y 1 - w , x y 0.5];ops sdpsettings(verbose, 3, solver, cplex);
sol optimize(Constraints,objective,ops);if sol.problem 0disp(求解成功);
elsedisp([求解失败错误原因为,sol.info]);
end 运行结果 显然这两个鲁棒优化无法使用Yalmip鲁棒优化模块求解的原因是一样的。都是任意一种最恶劣的场景都具有可行解但是无法在所有可能的恶劣场景下具有可行解。 了解了问题的本质之后我们可以知道Yalmip鲁棒优化模块的求解结果和我们的要求其实是不太一样的。但如果想使用Yalmip鲁棒优化模块求解需要做如下考虑 问题其实就是当仓库的总容量时肯定是可以满足任意一种恶劣的场景(即对于变量d任意的取值都可以得到可行解)但无法满足所有可能的场景(即使它们不会同时发生就和例子中的w0和w1不会同时发生一样)。如要考虑所有可能场景则需要仓库的总容量更大即 将这个约束条件到确定性优化中得到一个可行解其中z的取值为[546;0;274]再带入子问题中进一步得到求解结果 我们得到了鲁棒优化模块的求解结果但这个结果和文献[1]中所提的优化问题并不完全等价因此求解两阶段鲁棒优化时并不推荐使用这种方法。假设我们要做的是单阶段鲁棒优化使用Yalmip的鲁棒优化模块可以快速得到结果。例如将两阶段鲁棒优化改写成单阶段鲁棒优化 求解这个问题的matlab代码如下
%% 单阶段鲁棒优化问题
clc
clear
close all
warning off
yalmip(clear)%% 参数矩阵
f [400; 414; 326];
a [18; 25; 20];
k 800;
C [22, 33, 24;33, 23, 30;20, 25, 27];
d_ [206; 274; 220];
d_wave 40;
gamma [1.8,1.2];
P [1 1;1 1;1 0];%% 决策变量
y binvar(3,1);
z sdpvar(3,1);
x sdpvar(3,3,full);
g sdpvar(3,1);%% 目标函数
objective f*y a*z sum(sum(C.*x));%% 约束条件
Constraints [];
Constraints [Constraints , z 0 , x 0 ];
Constraints [Constraints , z k*y];
Constraints [Constraints , sum(x,2) z];
Constraints [Constraints , sum(x) (d_ g*d_wave)];G [uncertain(g) , g 0 , g 1 , g*P gamma];%% 设置求解器
opssdpsettings(verbose, 3, solver, gurobi);
soloptimize([Constraints , G],objective,ops);%% 分析错误标志
if sol.problem 0disp(求解成功);
elsedisp([求解失败错误原因为,sol.info]);
end%% 输出结果
if sol.problem 0objective value(objective)y value(y)x value(x)z value(z)
end 运行结果如下 另外由于使用KKT条件或强对偶变换的方式求解子问题后都要采用相同的方式和主问题迭代得到两阶段鲁棒优化问题的最优解所以我会先把两种求解方法介绍完之后再讲解如何将主问题和子问题结合起来迭代求解。
2.2 KKT条件求解子问题 为了方便求解我们首先把子问题的内层min优化问题写出来并将所有约束写成≤的形式 由于内层优化中相当于只有变量x不确定变量g和第一阶段优化变量y可以视为常数其拉格朗日函数为 其中该优化问题的最后三项为互补松弛条件是非线性约束可以引入0-1变量利用大M法进行等效线性化过程如下 其中q为m维0-1变量s为n维0-1变量v为m×n维0-1变量。将内层优化的KKT方程组添加到外层优化中就可以将双层优化问题转为单层优化问题如下式所示 经过上面的处理便顺利将max-min形式的子问题转为混合整数线性规划问题并可以使用Yalmip进行求解代码在压缩包中的Problem1文件夹中运行Problem1_subproblem_KKT.m文件即可得到结果(假设z[772;0;0])运行结果如下 使用KKT条件求解和使用鲁棒优化模块的结果略微有些不同同时使用KKT条件求解时可以返回不确定变量的取值即最恶劣场景下用户的需求分别为20627440314和2200.8*40252。
2.3 使用对偶变换求解子问题 子问题的内层优化问题为 其中αβγ均为对偶变量α为3×1的变量β为1×3的变量γ为3×3的变量我们以挨个条件对比的方式来写对偶问题 1)原问题为min问题对偶问题为max问题。 2)原问题中有9个变量因此对偶问题中有9个约束条件。 3)原问题中变量都≥0因此对偶问题中约束条件的符号是≤0。 4)原问题中目标函数的系数cij因此对偶问题的约束条件中≤号右边的常数也为cij。 5)原问题中有33915个约束条件因此对偶问题有15个决策变量。 6)原问题中约束条件的符号为≥因此对偶问题中决策变量的取值都≥0。 7)原问题中第1个约束条件中x11x12x13的系数均为1因此对偶变量α1在对偶问题9个不同的约束条件中的系数分别为[1,1,1;0,0,0;0,0,0]原问题中第2个约束条件中x21x22x23的系数均为-1因此对偶变量α2在对偶问题9个不同的约束条件中的系数分别为[0,0,0;1,1,1;0,0,0]原问题中第3个约束条件中x31x32x33的系数均为1因此对偶变量α3在对偶问题9个不同的约束条件中的系数分别为[0,0,0;0,0,0;1,1,1]。 8)原问题中第4个约束条件中x11x21x31的系数均为1因此对偶变量β1在对偶问题9个不同的约束条件中的系数分别为[1,0,0;1,0,0;1,0,0]原问题中第5个约束条件中x12x22x32的系数均为1因此对偶变量β2在对偶问题9个不同的约束条件中的系数分别为[0,1,0;0,1,0;0,1,0]原问题中第6个约束条件中x13x23x33的系数均为1因此对偶变量β3在对偶问题9个不同的约束条件中的系数分别为[0,0,1;0,0,1;0,0,1]。 9)原问题中第7-15个约束条件中xij的系数均为1因此对偶变量γij在对偶问题9个不同的约束条件中的系数均为1。 综上所述子问题内层优化的对偶问题可以写做 将内层优化的对偶问题和外层问题合并得到 其中目标函数中包含变量β和变量g的乘积是一个非线性项文献[1]的附录中假设变量g为0-1变量则可以使用大M法进行线性化。但这种假设需要满足一定的条件即不确定预算Γ只能为整数。当不确定预算Γ不是整数时只能将优化问题看作一个二次规划问题虽然是非线性优化但也可以使用KKT条件或者求解器求解这里我们直接使用求解器进行求解。 代码在压缩包中的Problem1文件夹中运行Problem1_subproblem_dual.m文件即可得到结果(假设z[772;0;0])运行结果如下 2.4 CCG算法KKT条件求解 2.1节到2.3节分别采用Yalmip工具箱的鲁棒优化模块KKT条件与对偶变换三种方法求解得到子问题。其中鲁棒优化模块求解与子问题的逻辑有一些区别。又因为子问题的对偶变换得到的是一个二次规划问题求解速度略慢于KKT条件得到的线性规划问题因此我们首先使用KKT条件求解子问题并采用CCG方法迭代求解两阶段鲁棒优化问题(如果想使用对偶变换步骤基本相同)。 首先可以把两阶段鲁棒优化问题可以分为主问题和子问题 KKT版本主问题MP1_KKT 网上很多CCG的代码其实对CCG算法的原理都没有理解透彻在每次迭代过程中都是在主问题中更新而不是增加决策变量和约束条件很容易出现迭代无法收敛的情况就算收敛了得到的也不是最优解。 KKT条件版本子问题SP1_KKT 使用CCG算法与KKT条件求解两阶段鲁棒优化的步骤概括如下 对比可知当子问题有解时需要将第k次迭代时子问题的目标函数作为约束条件添加到主问题中。如果无解时则需要相应地减少向主问题条件的约束条件。 使用KKT条件CCG算法求解该两阶段鲁棒优化问题的代码在压缩包中的Problem1文件夹中运行Problem1_KKT.m文件即可得到结果运行结果如下 结果表明该方法可以有效求解两阶段鲁棒优化问题。
2.5 CCG算法对偶变换求解 首先可以把两阶段鲁棒优化问题可以分为主问题和子问题 对偶变换版本主问题MP1_dual 对偶变换版本子问题SP1_dual 使用CCG算法与对偶变换求解两阶段鲁棒优化的步骤概括如下 使用对偶变换CCG算法求解该两阶段鲁棒优化问题的代码在压缩包中的Problem1文件夹中运行Problem1_dual.m文件即可得到结果运行结果如下 由于对偶变换引入了非线性项所以求解效率明显低于KKT条件的求解方式。
3.微电网两阶段鲁棒优化调度编程实战 更多内容请关注MatlabYalmip两阶段鲁棒优化通用编程指南(下)鲁棒优化入门(7)—MatlabYalmip两阶段鲁棒优化通用编程指南(下)
参考文献
[1]Zeng B, Zhao L. Solving two-stage robust optimization problems using a column-and-constraint generation method[J]. Operations Research Letters, 2013, 41(5): 457-461
[2]刘一欣,郭力,王成山.微电网两阶段鲁棒优化经济调度方法[J].中国电机工程学报,2018,38(14):4013-40224307.
PS 完整资料可以私信博主获取。
