0.引言
????????上一篇博客介紹了使用Yalmip工具箱求解單階段魯棒優(yōu)化的方法。這篇文章將和大家一起繼續(xù)研究如何使用Yalmip工具箱求解兩階段魯棒優(yōu)化(默認(rèn)看到這篇博客時(shí)已經(jīng)有一定的基礎(chǔ)了，如果沒有可以看看我專欄里的其他文章)。關(guān)于兩階段魯棒優(yōu)化與列與約束生成算法的原理，之前的博客已經(jīng)詳細(xì)地介紹過了，這里就不再過多介紹，主要是結(jié)合實(shí)例來講解編程思路。這篇博客用到了兩個(gè)算例，1個(gè)是兩階段魯棒優(yōu)化問題和列與約束生成算法的開山鼻祖[1]，另一個(gè)是電氣專業(yè)中兩階段魯棒優(yōu)化問題最熱門的文章之一[2]，相信大家在網(wǎng)上見到過無數(shù)號(hào)稱完美復(fù)現(xiàn)的代碼，但實(shí)際上大部分都是有問題的(包括我自己早期寫的代碼，也是被網(wǎng)上的代碼帶歪了，后面理解慢慢深入才發(fā)現(xiàn)問題所在)。
????????求解兩階段魯棒優(yōu)化問題一共有兩個(gè)難點(diǎn)，一是求解max-min或者min-max形式的子問題，其實(shí)就是求解一個(gè)單階段魯棒優(yōu)化，上一篇博客我已經(jīng)非常詳細(xì)地介紹了求解方式，借助Yalmip工具箱，共有三種不同的方式可以解決。二是主問題和子問題的迭代求解，也就是列與約束生成算法(C&CG)的實(shí)現(xiàn)。很多代碼在復(fù)現(xiàn)C&CG算法時(shí)并沒有向主問題同時(shí)添加列(變量)和約束，這也是代碼中最常見的問題。針對(duì)這兩個(gè)難點(diǎn)，我將用兩個(gè)不同的算例詳細(xì)地進(jìn)行講解。
????????此外，文獻(xiàn)[1]和[2]中都是采用了先將約束條件寫成緊湊的矩陣形式，然后再對(duì)子問題進(jìn)行處理的方式，很多朋友和我反映這部分太難處理了，實(shí)際問題的約束建模過程中經(jīng)常包括循環(huán)語句，想要轉(zhuǎn)成矩陣形式確實(shí)很不容易。這篇文章中我將分別采用兩種不同的方式求解魯棒優(yōu)化。一是采用原始的約束條件，省去將約束條件轉(zhuǎn)為矩陣形式的步驟，這種方式數(shù)學(xué)公式可能會(huì)更繁瑣，但比起矩陣形式的轉(zhuǎn)換，理解起來會(huì)更容易一些。二是采用矩陣形式進(jìn)行編程，在博客中我教大家一種非常簡單就能將約束寫為矩陣形式的方法，文中只是介紹了如何使用，之后也會(huì)單獨(dú)寫博客對(duì)此詳細(xì)展開。
????????總之，這篇博客干貨滿滿，可以認(rèn)真通讀一遍，跟著博客中的思路親自動(dòng)手使用Matlab+Yalmip實(shí)現(xiàn)兩階段魯棒優(yōu)化的編程(博客中提到的所有例子我都提供了相應(yīng)的代碼)。相信大家理解后，面對(duì)任何類型的兩階段魯棒優(yōu)化問題都能迅速使用類似的方法進(jìn)行解決。
????????博客中主要包含8大內(nèi)容：
????????①.拿到一個(gè)復(fù)雜的兩階段魯棒優(yōu)化問題的分析步驟和方法。
????????②.采用Yalmip工具箱中的uncertain函數(shù)和魯棒優(yōu)化模塊求解兩階段魯棒優(yōu)化的子問題。
????????③.Yalmip工具箱中的魯棒優(yōu)化模塊和常規(guī)的求解思路有什么異同。
????????④.使用KKT條件求解兩階段魯棒優(yōu)化的子問題，并使用C&CG算法進(jìn)行迭代求解。
????????⑤.使用對(duì)偶變換求解兩階段魯棒優(yōu)化的子問題，并使用C&CG算法進(jìn)行迭代求解。
????????⑥.采用Yalmip工具箱的內(nèi)置函數(shù)，將線性約束寫成緊湊矩陣形式的方法。
????????⑦.矩陣形式的兩階段魯棒優(yōu)化問題，如何快速寫出子問題內(nèi)層優(yōu)化的KKT條件，并使用C&CG算法進(jìn)行迭代求解。
????????⑧.矩陣形式的兩階段魯棒優(yōu)化問題，如何快速寫出子問題內(nèi)層優(yōu)化的對(duì)偶問題，并使用C&CG算法進(jìn)行迭代求解。
????????由于博客篇幅較長，將分上下兩篇發(fā)布，其中上篇使用的是文獻(xiàn)[1]中的算例，包含上述①-⑤的內(nèi)容。下篇使用文獻(xiàn)[2]中的算例，包含上述①、④、⑥-⑧的內(nèi)容。
????????這篇博客是下篇的內(nèi)容。
1.兩階段魯棒優(yōu)化基本形式
????????如文獻(xiàn)[1]中所述，標(biāo)準(zhǔn)的兩階段魯棒優(yōu)化問題的形式為：
????????其中，y為第一階段決策變量，u為不確定變量，x為第二階段決策變量。和分析單階段魯棒優(yōu)化問題的五個(gè)特征一樣，拿到一個(gè)復(fù)雜的兩階段魯棒優(yōu)化問題先不用慌，按照下面的步驟進(jìn)行分析即可：
????????1）確定第一階段決策變量有哪些，將其與變量y對(duì)應(yīng)。
????????2）確定第二階段決策變量有哪些，將其與變量x對(duì)應(yīng)。
????????3）確定不確定變量有哪些，將其與變量u對(duì)應(yīng)。
????????4）確定優(yōu)化問題中不確定集合的形式，并考慮是否可以直接使用Yalmip中的魯棒優(yōu)化模塊進(jìn)行求解。
????????5）確定目標(biāo)函數(shù)是否有僅包含第一階段決策變量的項(xiàng)，如果有的話可以單獨(dú)拿出來。
????????6）確定子問題的目標(biāo)函數(shù)，將其與魯棒優(yōu)化的標(biāo)準(zhǔn)形式相對(duì)應(yīng)。
????????7）確定約束條件，考慮是否包含非線性約束，是否需要線性化。
????????8）求解max-min或者min-max類型的子問題。
????????9）使用迭代方式，將子問題產(chǎn)生的變量和約束不斷添加到主問題中，最終得到最優(yōu)解。
????????下面分別以文獻(xiàn)[1]和[2]中的優(yōu)化問題進(jìn)行講解說明：
2.兩階段魯棒運(yùn)輸問題編程實(shí)戰(zhàn)
????????更多內(nèi)容，請(qǐng)關(guān)注Matlab+Yalmip兩階段魯棒優(yōu)化通用編程指南(上)：
魯棒優(yōu)化入門(6)—Matlab+Yalmip兩階段魯棒優(yōu)化通用編程指南(上)
3.微電網(wǎng)兩階段魯棒優(yōu)化調(diào)度編程實(shí)戰(zhàn)
????????第二個(gè)算例來源于文獻(xiàn)[2]，其兩階段魯棒優(yōu)化問題的形式為：
s.t.
????????我之前寫過一篇博客解析這篇論文，但是使用的是緊湊的矩陣形式。這次博客我先嘗試使用一般形式的約束進(jìn)行求解，方便大家體會(huì)采用矩陣形式求解兩階段魯棒優(yōu)化的方便之處。
????????首先逐步進(jìn)行分析：
????????1）確定第一階段決策變量有哪些。
????????2）確定第二階段決策變量有哪些。
????????3）確定不確定變量有哪些。
????????4）確定優(yōu)化問題中不確定集合的形式，并考慮是否可以直接使用Yalmip中的魯棒優(yōu)化模塊進(jìn)行求解。
????????該優(yōu)化問題的不確定集中為箱式不確定集，但包含0-1變量，不可以直接使用Yalmip中的魯棒優(yōu)化模塊進(jìn)行求解。
????????5）確定目標(biāo)函數(shù)是否有僅包含第一階段決策變量的項(xiàng)，如果有的話可以單獨(dú)拿出來。
????????該優(yōu)化問題第一階段和第二階段目標(biāo)函數(shù)相同。
????????6）確定子問題的目標(biāo)函數(shù)，將其與兩階段魯棒優(yōu)化的標(biāo)準(zhǔn)形式相對(duì)應(yīng)。
????????7）確定約束條件，考慮是否包含非線性約束，是否需要線性化。
????????子問題中的約束條件均為線性，但不確定變量的約束條件包含0-1變量，則含變量的優(yōu)化問題并不滿足強(qiáng)對(duì)偶定理，KKT條件和對(duì)偶變換都無法適用。那是不是無法使用常規(guī)的方法求解單階段魯棒優(yōu)化形式的子問題呢？
????????當(dāng)然不是，如果出現(xiàn)這種情況，我們就可以把不確定變量寫在外層。那么對(duì)于子問題的內(nèi)層優(yōu)化問題來說，變量是一個(gè)確定值，即使其中包含0-1變量，也只是一個(gè)常數(shù)。實(shí)際上子問題的內(nèi)層優(yōu)化并不包含0-1變量，因此還是可以使用常規(guī)的KKT條件或強(qiáng)對(duì)偶定理進(jìn)行求解。但是由于不確定變量中包含0-1變量，直接通過uncertain函數(shù)直接使用Yalmip魯棒優(yōu)化模塊的方法不可行。
????????分析完成后，分別使用KKT條件和強(qiáng)對(duì)偶定理求解子問題，具體如下：
3.1?KKT條件求解子問題
????????為了方便求解，我們首先把子問題的內(nèi)層min優(yōu)化問題寫出來，并將所有約束寫成≤0的形式或=0的形式：
????????根據(jù)拉格朗日函數(shù)可以寫出內(nèi)層優(yōu)化的KKT條件。
????????其中表示上三角全為1，主對(duì)角線以下全為0的24階方陣，即：
????????要想求解這個(gè)問題，還需要寫出引入16組0-1變量，使用大M法將16組互補(bǔ)松弛條件轉(zhuǎn)為線性約束，具體如下：
????????其中q1-q16都是0-1變量，qi和λi的維度相同。
????????將內(nèi)層優(yōu)化的KKT方程組添加到外層優(yōu)化中，就可以將雙層優(yōu)化問題轉(zhuǎn)為單層優(yōu)化問題。經(jīng)過上面的處理，便順利將max-min形式的子問題轉(zhuǎn)為混合整數(shù)線性規(guī)劃問題，并可以使用Yalmip進(jìn)行求解，代碼在壓縮包中的Problem2文件夾中，運(yùn)行Problem2_subproblem_KKT1.m文件即可得到結(jié)果。
????????PS：假設(shè)第一階段變量的取值：
????????運(yùn)行結(jié)果如下：
????????從運(yùn)行結(jié)果可以看出，使用KKT條件可以正常求出最優(yōu)解，但由于引入了非常多的拉格朗日乘子和相應(yīng)的0-1變量，對(duì)于初始的約束形式手動(dòng)寫KKT條件非常麻煩，一旦出錯(cuò)也很難發(fā)現(xiàn)問題在哪。
3.2 將兩階段魯棒優(yōu)化問題寫成矩陣形式
????????對(duì)于約束條件比較復(fù)雜的兩階段魯棒優(yōu)化問題，先把目標(biāo)函數(shù)和約束條件寫成緊湊的矩陣形式再使用KKT條件或?qū)ε甲儞Q求解會(huì)方便很多。文獻(xiàn)[1]中將優(yōu)化問題寫成了緊湊的矩陣形式，但存在一定問題，我在此重新梳理一遍。
????????首先，將優(yōu)化問題寫成緊湊的矩陣形式(等式約束都改寫成不等式形式，例如x=0可以寫成x≥0，-x≥0)：
????????其中，第一階段決策變量x是一個(gè)48維的列向量：
????????第二階段決策變量y是一個(gè)192維的列向量：
????????c和y的維度一樣，但c是一個(gè)192維的常數(shù)列向量：
????????那么如何確定矩陣G，矩陣E，矩陣M，向量h的取值呢？之前的博客中我采用了手動(dòng)推導(dǎo)的方式(微電網(wǎng)兩階段魯棒優(yōu)化經(jīng)濟(jì)調(diào)度方法matlab代碼_)。
????????這次我來教大家使用Yalmip工具箱中的函數(shù)depends、getbase、getbasematrix、see寫出約束矩陣取值的方法，函數(shù)的語法可以參考官方文檔。
????????下面是一個(gè)簡單的線性規(guī)劃問題：
????????我用這個(gè)例子來幫助大家體會(huì)depends、getbase、getbasematrix、see等函數(shù)的用法，代碼如下：
clc
clear
close all
warning off
yalmip('clear')
%% 決策變量
sdpvar x1 x2
%% 目標(biāo)函數(shù)
obj = x1 + 2*x2;
z = [-2*x1 + 3*x2 - 12 ;
x1 + x2 - 14;
-3*x1 + x2 + 3;
3*x1 + x2 - 30];
x1_index = depends(x1)
x2_index = depends(x2)
M1 = full(getbase(z))
M2 = full(getbasematrix(z,depends(x1)))
M3 = full(getbasematrix(z,depends(x2)))
see(z)????????運(yùn)行結(jié)果如下：
????????觀察可知，depends函數(shù)返回的是變量在Yalmip工具箱中的編號(hào)，getbase函數(shù)以稀疏矩陣的形式返回sdpvar變量中的常數(shù)項(xiàng)系數(shù)以及變量系數(shù)，getbasematrix函數(shù)以稀疏矩陣的形式返回sdpvar變量中指定編號(hào)的變量系數(shù)，see函數(shù)直接在命令行輸出sdpvar變量中的常數(shù)項(xiàng)、變量系數(shù)以及用到的變量編號(hào)。
????????假設(shè)我們想把上面的優(yōu)化問題寫成緊湊的矩陣形式：
????????求矩陣A、b、c，然后使用矩陣形式求解優(yōu)化問題的代碼如下：
clc
clear
close all
warning off
yalmip('clear')
%% 決策變量
sdpvar x1 x2
%% 求系數(shù)矩陣
obj0 = x1 + 2*x2;
z = [-2*x1 + 3*x2 - 12 ;
x1 + x2 - 14;
-3*x1 + x2 + 3;
3*x1 + x2 - 30];
M1 = full(getbase(z));
M2 = full(getbase(obj0));
index = depends([x1 x2]);
A = M1(:,index + 1);
b = -M1(:,1);
c = M2(index + 1)';
%% 矩陣形式的目標(biāo)函數(shù)
x = [x1;x2];
obj = c'*x;
C = [A*x <= b , x >= 0];
%% 求解優(yōu)化問題
ops = sdpsettings('verbose', 3, 'solver', 'gurobi');
sol = optimize(C , -obj ,ops);
%% 判斷求解是否成功
if sol.problem == 0
disp('求解成功!!!');
x = value(x)
else
disp(['求解失敗，原因?yàn)?,sol.info]);
end????????運(yùn)行結(jié)果如下：
????????采用相同的方法，將文獻(xiàn)[1]所提確定性優(yōu)化問題改寫成矩陣形式求解的代碼在壓縮包中的Problem2文件夾中，運(yùn)行Problem2_matrix.m文件即可得到結(jié)果。
3.3 矩陣形式的C&CG算法與KKT條件求解
????????原始子問題為：
????????內(nèi)層優(yōu)化中變量x和不確定變量u都可以看作常數(shù)，其拉格朗日函數(shù)可以寫作：
????????其中，π和θ都是拉格朗日乘子。根據(jù)拉格朗日函數(shù)進(jìn)一步寫出KKT條件：
????????這就是初始的KKT條件，可以進(jìn)一步化簡。首先根據(jù)式(1)和(4)可以得到：
????????然后，根據(jù)式(1)和(3)可以得到：
????????經(jīng)過這樣處理，就可以消去變量θ，減少?zèng)Q策變量的數(shù)目，得到：
????????其中包含兩項(xiàng)互補(bǔ)松弛條件，可以通過引入0-1變量，使用大M法轉(zhuǎn)換為線性約束：
????????將KKT條件添加到子問題的外層優(yōu)化中，子問題便轉(zhuǎn)為：
????????子問題變成了一個(gè)混合整數(shù)線性規(guī)劃，可以采用求解器有效地進(jìn)行求解。
????????再復(fù)習(xí)一下，我們把文獻(xiàn)[2]中所提兩階段魯棒優(yōu)化分解成一個(gè)主問題MP2_KKT和一個(gè)子問題SP2_KKT：
????????主問題MP2_KKT：
????????子問題SP2_KKT：
????????使用KKT條件+C&CG算法求解該兩階段魯棒優(yōu)化問題的代碼在壓縮包中的Problem2文件夾中，運(yùn)行Problem2_KKT.m文件即可得到結(jié)果。
3.4 矩陣形式的C&CG算法與對(duì)偶變換求解
????????原始子問題為：
????????其內(nèi)層優(yōu)化min問題的對(duì)偶問題可以寫做：
????????將其與外層max問題合并得到：??
????????該問題的目標(biāo)函數(shù)中包含的非線性項(xiàng)，可以進(jìn)行線性化，也可以直接使用求解器進(jìn)行求解。為了方便起見，我直接使用了求解器進(jìn)行求解。但是由于問題規(guī)模比較大，求解時(shí)間會(huì)比較長(大概要3小時(shí)左右才能得到最優(yōu)解)。
????????首先，將不確定變量u寫成：
????????式中，符號(hào)表示矩陣的Hadamard乘積(哈達(dá)瑪積 Hadamard Product - 知乎)，即對(duì)應(yīng)元素相乘，對(duì)應(yīng)于matlab中的“.*”符號(hào)。在此基礎(chǔ)上，可以對(duì)原優(yōu)化問題的目標(biāo)函數(shù)進(jìn)行化簡：
????????兩階段魯棒優(yōu)化問題可以分為主問題和子問題：
????????對(duì)偶變換版本主問題MP2_dual：
????????對(duì)偶變換版本子問題SP2_dual：
????????使用C&CG算法與對(duì)偶變換求解兩階段魯棒優(yōu)化的步驟概括如下：
????????使用對(duì)偶變換+C&CG算法求解該兩階段魯棒優(yōu)化問題的代碼在壓縮包中的Problem2文件夾中，運(yùn)行Problem2_dual.m文件即可得到結(jié)果。
參考文獻(xiàn)：
[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]劉一欣,郭力,王成山.微電網(wǎng)兩階段魯棒優(yōu)化經(jīng)濟(jì)調(diào)度方法[J].中國電機(jī)工程學(xué)報(bào),2018,38(14):4013-4022+4307.文章來源：http://www.zghlxwxcb.cn/news/detail-793736.html
PS：
????????完整資料可以私信博主獲取。文章來源地址http://www.zghlxwxcb.cn/news/detail-793736.html
到了這里，關(guān)于魯棒優(yōu)化入門(7)—Matlab+Yalmip兩階段魯棒優(yōu)化通用編程指南(下)的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
