SeDuMi詳細(xì)教程

1、SeDuMi User Guide線性規(guī)劃（LP問題：Linear Programming）將需要解決的問題表述成原始問題：或者對偶問題：例子：為了解決該LP問題，我們必須添加一些松弛變量，比如x3和x4，我們即可在MATLAB 中輸入b和c向量以及A矩陣:現(xiàn)在我們就可以通過sedumi指令來求解這個問題 sedumi(A,b,c)eqs m = 2, order n = 5, dim = 5, blocks = 1 nnz(A) = 6 + 0, nnz(ADA) = 4, nnz(L) = 3it :b*ygapdelta rate t/tP* t/tD* feas cg cg prec0

2、 : 6.02E+000 0.0001 : -1.20E+000 1.84E+000 0.000 0.3055 0.9000 0.9000 1.86 1 12.1E+0002 : -6.94E-002 5.15E-001 0.000 0.2803 0.9000 0.9000 2.14 1 1 8.6E-0013 : -1.21E-001 2.72E-002 0.000 0.0528 0.9900 0.9900 1.16 1 1 4.2E-0024 : -1.25E-001 8.69E-006 0.000 0.0003 0.9999 0.9999 1.02 1 1iter seconds dig

3、itsc*xb*y40.3Inf-1.2500000000e-001 -1.2500000000e-001|Ax-b| = 0.0e+000, Ay-c_+ = 1.7E-016, |x|= 2.9e+000, |y|= 2.8e-001Detailed timing (sec)PreIPMPost1.404E-001 2.964E-001 4.680E-002 Max-norms: |b|=5, |c| = 1, Cholesky |add|=0, |skip| = 0, |L.L| = 1. ans =(1,1)1.9583(2,1)2.0833表明最優(yōu)值為-1.2500000000e-0

4、01=-0.125，對應(yīng)c*x以下的值，同時返回最優(yōu)解：x1=1.9583,x2=2.0833,發(fā)現(xiàn)x確實有解，因為其每一個元素都是非負(fù)的，而且Ax=b，可以用命令min(x)和norm(Ax-b)來檢驗。當(dāng)然，也會出現(xiàn)一些舍入誤差，如下： norm(A*x-b) ans =1.7764e-015 norm(A*(24*x)-24*b)ans = 0二次型和半定限制（Quadratic and semi definite constraint）在sedumi中，有可能強(qiáng)制加入二次型或者半定限制條件，即通過限制變量進(jìn)入一個二次型和核或者半正定矩陣的核，這樣一個限制條件代替了線性規(guī)劃中的非負(fù)性條件

5、，需要x屬于K，一個對稱核中，它是一個非負(fù)象限,二次型核和半正定矩陣核的笛卡兒積。最優(yōu)化問題的標(biāo)準(zhǔn)形式為：對偶形式為：二次型核二次型核由以下形式來確定：考慮以下問題：其中P為給定矩陣，q為給定向量，以上是魯棒最小均方問題。其中決策變量為標(biāo)量y1和y2，還有向量y3.這個問題有兩個二次型限制：給定P和q，以下MATLAB函數(shù)（rls.m）將該問題表述成標(biāo)準(zhǔn)對偶問題。A矩陣被表述為轉(zhuǎn)置方向，即表述為At。Rls.m% At,b,c,KI = rls(P,q)% Creates dual standard form for robust least squares problem Pu=q.func

6、tion At ,b,c,K= rls(P,q) m, n = size(P) ;% - - - - - - - - - - minimize y-1 + y-2 -b = -sparse(1; 1; zeros(n,1);% - (y_1, q - p y_3) in Qcone -At = sparse (-1, zeros(1,1+n); .zeros(m, 2), P); c = 0;q ;K.q=1+m;%- (y_2, (1, y_3) in Qcone -At = At; 0, -1, zeros(1,n); .zeros(1,2+n); .zeros(n,2),-eye(n);

7、c = sparse(c;0;1;zeros(n,1);K.q= K.q, 2+n;我們可以注意到以上的方程中運(yùn)用的是系數(shù)數(shù)據(jù)存儲形式，為的是節(jié)省內(nèi)存。此外，還定義了一個結(jié)構(gòu)體K，其中K.q列出了二次型核的維數(shù)。K結(jié)構(gòu)體將被運(yùn)用于告知SeDuMi：c-ATy的元素并不被限制為非負(fù)當(dāng)他們都被用于線性規(guī)劃中。反而言之，第一個行K.q(1)被限制在一個二次型核中，最后一行K.q(2)被限制在另外一個二次型核中，這就是我們在之前建立對稱核的方法，而且建立了兩個二次型限制。作為數(shù)值仿真的例子，我們來來解決一個魯棒最小平方值的問題，其依靠P中的列 P = 3 1 4; 0 1 1;-2 5 3; 1 4

8、5; q = 0; 2; 1; 3; At,b,c,K = rls(P,q); x,y,info = sedumi(At,b,c,K);運(yùn)行后出現(xiàn)的結(jié)果是SeDuMi 1.1R3 by AdvOL, 2006 and Jos F. Sturm, 1998-2003.Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta =0.500eqs m = 5, order n = 5, dim = 11, blocks = 3 nnz(A) = 16 + 0, nnz(ADA) = 23, nnz(L) = 1

9、4it :b*ygapdelta ratet/tP* t/tD*feas cg cg prec0 :4.00E+000 0.0001 : -3.10E+000 3.01E-001 0.000 0.0751 0.9900 0.99001.55 1 1 3.4E-0012 : -3.33E+000 6.02E-003 0.000 0.0200 0.9900 0.99001.03 1 1 6.2E-0033 : -3.33E+000 1.31E-005 0.000 0.0022 0.9990 0.99901.00 1 1 1.3E-0054 : -3.33E+000 1.26E-006 0.367

10、0.0958 0.9900 0.99001.00 1 1 1.3E-0065 : -3.33E+000 3.81E-008 0.000 0.0303 0.9900 0.99001.00 1 1 3.8E-0086 : -3.33E+000 2.24E-009 0.425 0.0587 0.9902 0.99001.00 1 1 2.3E-009iter seconds digitsc*xb*y6 0.0 9.4 -3.3329085951e+000 -3.3329085963e+000|Ax-b| = 9.5e-010, Ay-c_+ = 2.1E-009, |x|= 2.0e+000, |y

11、|= 2.5e+000Detailed timing (sec)PreIPMPost3.120E-0023.120E-0020.000E+000Max-norms: |b|=1, |c| = 3,Cholesky |add|=0, |skip| = 0, |L.L| = 1.在以上這些SeDuMi的語句中，我們發(fā)現(xiàn)了一個新的輸入量viz.K，這個量使得SeDuMi能夠解決一個（5）和（6）形式的最優(yōu)化問題，其中對稱核K被描述為結(jié)構(gòu)體K。沒有K的話，SeDuMi只能解決（1）和（2）的問題。我們可以通過驗證（7）中的不等式來確定結(jié)果y是否能夠滿足條件（9）。然而，SeDuMi提供了一種更簡單的方

12、法eigK，這個函數(shù)返回相對應(yīng)于對稱核的矩陣的特征值。一個對稱形核包括哪些僅有非負(fù)特征值的向量，其中參見Faraut and Koranyi9。如下，我們能因此檢驗可行性和最優(yōu)性： eigK(x, K), eigK(c-At*y, K) ans =0.0000 -0.00001.4142 3.23070.0000 -0.00001.4142 1.4827 x*(c-At*y)ans =3.4744e-009對于對稱核K來說，總是有xTz=0，對所有屬于K的z和x成立，因此，x提供里一種在線性規(guī)劃中對y最優(yōu)化的認(rèn)證。Farkas的對偶方案的分解也是采用同樣的思想。具體參見15的方法。然而，一個奇

13、怪的現(xiàn)象出現(xiàn)了，x和y幾乎是可解得，然而cTx-bTy大多是負(fù)的（|x|和或者|y|然后肯定會是非常大的），根據(jù)原理，SeDuMi將會返回一個無窮大的精確的digits。這個現(xiàn)象在14和23中也有解釋。我們需要讓這個最優(yōu)化模型有著非負(fù)的和二次型的核限制，而這一點(diǎn)是可行的。比如，我們可以以上例子中的限制y310來解決模型注意到該問題是無解的：1/x1下確界是0，對x1趨近與無窮大來說：clear K;c= 0,1,0;b = sqrt(2); A = 0, 0, 1; K.r = 3;x,y,info=sedumi(A,b,c,K);x(2), x(1)*x(2)ans =6.5278e-005

14、ans =1.0695你可能會發(fā)現(xiàn)x2根本沒有足夠趨近于0，而且x1也不是等于無窮大。然而，這個原始解是可解的，對偶解幾乎是可解的，而且對偶gap幾乎是負(fù)的。這就解釋了一個誤差范圍困難，而且對這種不規(guī)則的問題是很常見的。在Section 5 中，我們可以看看到底如何獲得一個更精確的解決方案，通過設(shè)計一個選項參數(shù)，pars.eps。上述例子可以解釋，K.r來列出Rcone限制的維數(shù)，類似于Qcone限制中K.q的定義，設(shè)定了K.l，K.q，K.r也就引出了以下形式的對稱核。舉例說明，我們可以添加界x1 X=vec(1,5,-3;5,2,-9;-3,-9,4) X =1 5 -3 5 2 -9 -

15、3 -9 4Vec的逆運(yùn)算就是mat。因此，如果x是一個n2長度的向量，而mat（x）就創(chuàng)建了一個nn的矩陣，然后依次填入x向量的元素，比如： mat(X) ans =1 5 -35 2 -9-3-94以下MATLAB函數(shù)產(chǎn)生了一個問題（11）的標(biāo)準(zhǔn)的原始形式：% At,b,c,K =specfac(b)% Creates primal standard form for minimal phase spectral factorization. function At,b,c,K =specfac(b)m = length(b) ;% - minimize sum (m-i)*X(i,i)

16、-c = vec(spdiags(m-1:-l:0),0,m,m);% - Let e be all-1, and allocate space for the A-matrix -e = ones(m,1);At = sparse( , , ,m2,m,m*(m+l)/2) ;% - sum(diag(): k) = b(k) -for k= 1:mAt(:,k) = vec(spdiags(e,k-l,m,m);endK.s = m;字段K.s=m告訴SeDuMi我們需要一個mm的矩陣mat（x），而且滿足對稱正定的。我們能解決如下的問題：b = 2; 0.2; -0.3;At, b, c

17、 ,K = specfac(b) ;x,y,info = sedumi(At,b,c,K);得到的結(jié)果為：it :b*ygapdelta ratet/tP* t/tD* feas cg cgprec0:1.69E+001 0.0001: -2.44E-001 3.91E+000 0.000 0.2312 0.9000 0.90001.39 1 12.0E+0002: 1.86E-001 3.22E-001 0.000 0.0824 0.9900 0.99001.31114.0E-0013: 1.23E-001 1.12E-004 0.000 0.0003 0.9999 0.99990.9911

18、2.6E-0044: 1.23E-001 5.18E-006 0.000 0.0463 0.9900 0.99001.00111.2E-0055: 1.23E-001 2.36E-007 0.000 0.0455 0.9900 0.99001.00116.8E-0076: 1.23E-001 1.90E-008 0.327 0.0806 0.9900 0.99001.00225.5E-0087: 1.23E-001 8.16E-010 0.000 0.0429 0.9906 0.99000.99222.2E-009iter seconds digitsc*xb*y7 0.2 8.5 1.227

19、3256559e-001 1.2273256524e-001|Ax-b| = 1.1e-010, Ay-c_+ = 1.4E-010, |x|= 2.0e+000, |y|= 7.6e-001Detailed timing (sec)PreIPMPost3.120E-002 2.028E-001 3.120E-002 Max-norms: |b|=2, |c| = 2,Cholesky |add|=0, |skip| = 0, |L.L| = 1.為檢驗正定性，可以采用MATLAB中的eig或者SeDuMi中的eigK： eig(mat(x),eigK(x,K) ans =0.0000 0.00000.0000 0.00002.00002.0000其中eigK更加方便，特別是當(dāng)有多個正定條件限制的時候，或者也有非負(fù)和二次型核限制。SeDuMI將總是產(chǎn)生對稱決策變量，比