1.等式約束二次規(guī)劃問題

(1)

其中??且為對(duì)稱矩陣，,?，. 求解等式約束二次規(guī)劃問題一般有直接消元法，正交分解法和拉格朗日法。

拉格朗日法

構(gòu)建Eq.(1)的拉格朗日乘子函數(shù)，得

由KKT條件可得：

寫成矩陣的形式得

系數(shù)矩陣稱為Lagrange矩陣，它是對(duì)稱不定的。如果Lagrange矩陣的逆存在

可得最優(yōu)解為.如果存在，可得

如果Lagrange矩陣的逆不存在，可以求偽逆。

2.等式和不等式約束二次規(guī)劃問題

由于帶有等式約束的二次規(guī)劃問題已經(jīng)有成熟的解法，所以要求解帶有不等式約束的二次規(guī)劃問題，一種很自然的想法就是將不等式約束轉(zhuǎn)換成等式約束。有效集法（Active Set Method）是求解帶不等式約束的二次規(guī)劃問題的一種經(jīng)典算法。首先看下問題的幾何表述

• 可行域：目標(biāo)問題的可行域是空間中的一個(gè)凸多面體（包括邊界和內(nèi)部）。例如上圖中的綠色多邊形區(qū)域。
• 等值面：如果矩陣G正定，那么目標(biāo)函數(shù)的等值面是空間中的一族同心超橢球面，所有橢球面相似且同向（各軸同向重合）。例如上圖中的同心橢圓曲線族。越內(nèi)層的橢球面，目標(biāo)函數(shù)在其上的取值越優(yōu)（低）。
• 優(yōu)化問題：在幾何表述下，我們的目標(biāo)就是在可行域內(nèi)找到達(dá)到同心橢圓曲線族最內(nèi)層的點(diǎn)。例如上圖中，P是無約束最優(yōu)解，Q是約束最優(yōu)解。

顯然加入約束之后，最優(yōu)解必定出現(xiàn)在某個(gè)或某幾個(gè)不等式約束的邊界上。

• 有效集：任給一個(gè)可行解，其有效集就是滿足等號(hào)成立的約束條件的集合。顯然，有效集必然包含所有等式約束，同時(shí)包含不等式約束的一個(gè)子集。注意這里是等號(hào)成立，也就是說，當(dāng)前解讓一個(gè)不等式約束的等號(hào)成立時(shí)，這個(gè)不等式約束才會(huì)被納入有效集。
• 最優(yōu)有效集：最優(yōu)解的有效集稱為最優(yōu)有效集。只需把最優(yōu)有效集中的約束全部改寫為等式約束，并扔掉其它約束，然后用Lagrange乘子法直接求解即可。

由于約束條件數(shù)目是有限的，而最優(yōu)有效集是其子集，進(jìn)而最優(yōu)有效集只有有限多種可能。于是很自然地就想到了下面大名鼎鼎的暴力破解法——窮舉法。

每次，我們從全體約束中選出一個(gè)子集（包含全部等式約束和部分不等式約束），稱之為本次嘗試的工作集，并求解該工作集對(duì)應(yīng)的子優(yōu)化問題。一個(gè)工作集對(duì)應(yīng)的子優(yōu)化問題是這樣定義的：在原問題中，將工作集中的約束全部改為等式約束，同時(shí)扔掉工作集之外的約束。由于這是一個(gè)等式約束問題，所以可以用Lagrange法輕松求解。之后檢查該解是否是原問題的可行解（即它是否也滿足工作集之外的約束）。如果是，則記錄下來，否則就扔掉。遍歷全部可能的子集之后，在記錄下來的所有解中，選一個(gè)使目標(biāo)函數(shù)值最優(yōu)的即可。

窮舉法雖然已經(jīng)可以確保在有限步內(nèi)獲得目標(biāo)問題的解，但其運(yùn)算量常常超乎想象（尤其是控制標(biāo)量和不等式約束條件個(gè)數(shù)較多時(shí)），所以我們?cè)谝韵聝煞矫鎸で蟾倪M(jìn)：

• 改進(jìn)最優(yōu)解的識(shí)別方式：利用凸二次規(guī)劃的特點(diǎn)，建立一組識(shí)別規(guī)則（實(shí)際上就是KKT條件），可直接判斷一個(gè)可行解是否是最優(yōu)解，因此算法一旦試出最優(yōu)解就立即停止，不必遍歷剩余可能性
• 改進(jìn)嘗試各種可能性的順序：不再隨機(jī)選取下一次嘗試的工作集，而是通過求解子優(yōu)化問題找到能使目標(biāo)函數(shù)值有效下降的新工作集和迭代點(diǎn)。

經(jīng)過以上兩項(xiàng)改進(jìn)，就形成了效率大幅提高的有效集法。

有效集法

初始可行解：有效集法需要一個(gè)初始可行解??作為迭代的起點(diǎn)。由于可行解與目標(biāo)函數(shù)無關(guān)，所以獲得初始可行解的方法與線性規(guī)劃完全一致，典型方法有兩階段法和大M法，，因此我們假設(shè)已經(jīng)獲得一個(gè)初始可行解 x0 ，并得到有效集為W0

迭代求解。迭代的過程主要包括兩個(gè)部分，一個(gè)尋找當(dāng)前有效集中的最優(yōu)解，另一個(gè)是更新有效集以得到最優(yōu)集。假設(shè)迭代k（k=1,2,3,...）步后，當(dāng)前解為xk，工作集為wk,根據(jù)xk落在可行域上的位置，可以劃分出四種情形。

• 情形1：?若xk??滿足KKT條件，則它就是最優(yōu)解，算法停止。
• 情形2：若 xk??在當(dāng)前工作集（等式約束集）下已達(dá)最優(yōu)，但工作集中存在至少一個(gè)約束，它原來是不等式約束，而且當(dāng) xk??朝該不等式約束的可行一側(cè)移動(dòng)時(shí)，可以使目標(biāo)函數(shù)值進(jìn)一步下降，那么我們令迭代點(diǎn)不動(dòng)（即 xk+1=xk并從工作集中去掉那個(gè)約束（即將它還原放松成不等式約束）。
• 情形3：若 xk在當(dāng)前工作集（等式約束集）下還可以改進(jìn)，但步長不能走滿（否則會(huì)走出原始問題的可行域），那么就取 xk?延改進(jìn)方向恰走到原始可行域邊界處的位置為下一個(gè)迭代點(diǎn) xk+1?，同時(shí)把碰到的邊界處對(duì)應(yīng)的不等式約束加入工作集。情形3的例子請(qǐng)參見下圖中x3?
• 情形4：若 xk在當(dāng)前工作集（等式約束集）下還可以改進(jìn)，并且步長能夠走滿（即改進(jìn)到最優(yōu)仍未出原始可行域），那么就取這個(gè)最優(yōu)點(diǎn)作為下一個(gè)迭代點(diǎn) xk+1，并保持工作集不變。情形4的例子請(qǐng)參見下圖中 x1和 x4。

下面是使用有效集求解帶不等式約束優(yōu)化問題的經(jīng)典例子。

G = eye(2,2)*2; d = [-2 -5]'; A = [-1 2; 1 2; 1 -2; -1 0; 0 -1]; b = [2;6;2;0;0]; Aeq = [-1 1]; beq = [0.03]; MaxIterations = 100; Tolerance = 1.0e-6; x0 = [0;1]; workset0= [3;5]; x = active_set_method(x0,G,d, A, b, Aeq, beq, MaxIterations, Tolerance) % options = optimoptions(@quadprog,'Algorithm','active-set','MaxIterations',500); y = quadprog(G, d, A, b, Aeq, beq, [], [], x0, options)function x = active_set_method(x0,G,d, Ai, bi, Aeq, beq, MaxIterations, Tolerance) %min 1/2*x'Gx+d'x %st. A*x <= b % Aeq*x = bne = size(Aeq,1);ni = size(Ai,1);workset = zeros(ni,1);x = x0;%給定初始解for i = 1:ni%將滿足初始解的不等式都加入有效集if abs(Ai(i,:)*x0 -bi(i)) <= Toleranceworkset(i) = 1;endendfor k = 1:MaxIterations%迭代求解 a = [Aeq; Ai(workset == 1,:)];b_ = [beq; bi(workset == 1,:)];[x1, lambda]=kkt(G, d, a, b_);s = x1 - x;if norm(s,2) < Tolerance%得到的解是本次有效集中的最優(yōu)解x = x1; min_lambda = 0;if length(lambda) > ne[min_lambda,index] = min(lambda(ne+1:end));endif min_lambda >= 0 %滿足KKT條件，得到最優(yōu)解break;else%尋找新的有效集for i=1:niif workset(i)&&sum(workset(1:i))==indexworkset(i)=0;break;endend endelse%得到的解不是本次有效集中的最優(yōu)解，繼續(xù)優(yōu)化[alpha, index] = Alpha(workset, x, s, Ai, bi);x = x + alpha*s;if alpha <1workset(index) = 1;endendendfunction [s, lambda]=kkt(G, g, Ae, be)ws = size(Ae,1);kkt_A = [G, Ae'; Ae, zeros(ws,ws)];kkt_B = [-g;be];slambda = pinv(kkt_A)*kkt_B; %kkt_A不一定可逆，需要求偽逆dim = size(G,2);s = slambda(1:dim);lambda = slambda(dim+1:end);endfunction [alpha, index] = Alpha(workset, x, s, Ai, bi)outset = find(workset == 0);Atp = Ai(outset,:)*s;Atp_index = Atp>0;out_plus_ind = outset(Atp_index);alphatset = (bi(out_plus_ind) - Ai(out_plus_ind,:)*x)./Atp(Atp_index);[alpha, index] = min([alphatset;1]);if alpha < 1index = out_plus_ind(index);endend end

參考：

https://zhuanlan.zhihu.com/p/29525367

https://zhuanlan.zhihu.com/p/26514613

https://blog.csdn.net/dymodi/article/details/50358278

https://zhuanlan.zhihu.com/p/50906541

總結(jié)

以上是生活随笔為你收集整理的二次规划--积极集法（active set method）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。