SeDuMi詳細(xì)教程

1、SeDuMi User Guide線性規(guī)劃(LP 問題：Lin ear Programmi ng )將需要解決的問題表述成原始問題:minimi2e rTJrsuch that Ax = bz, 0 fnr t - 113t. K K ,或者對(duì)偶問題:mxiniizesuch that G - 0 for i = 1 f 2,. 1 n例子:minimize T| -such that 10 xi 7x2 2 5xi + rj/2 0, r； 0.為了解決該LP問題，我們必須添加一些松弛變量，比如x3和x4,我們即可在MATLAB中輸入b和c向量以及A矩陣: c 二 ti； -1； o； o；

2、 A = 10. -了了. -1, 0, lt 1/21 0, 1; b - B; 3;現(xiàn)在我們就可以通過sedumi指令來求解這個(gè)問題sedumi(A,b,c)eqs m = 2, order n = 5, dim = 5, blocks = 1nn z(A) = 6 + 0, nn z(ADA) = 4, nn z(L) = 3feas1.7E-016, |x|= 2.9e+000, |y|=Post4.680E-002it :b*y gap delta rate t/tP* t/tD*cg cgprec0 :6.02E+000 0.0001 : -1.20E+000 1.84E+000

3、0.000 0.3055 0.9000 0.9000 1.86 1 1 2.1E+0002 : -6.94E-002 5.15E-001 0.000 0.2803 0.9000 0.90002.14118.6E-0013 : -1.21E-001 2.72E-002 0.000 0.0528 0.9900 0.99001.16114.2E-0024 : -1.25E-001 8.69E-006 0.000 0.0003 0.9999 0.99991.0211itersecondsdigits c*x b*y4 0.3 Inf-1.2500000000e-001-1.2500000000e-00

4、1|Ax-b| =0.0e+000, Ay-c_+ =2.8e-001Detailed timing (sec)Pre IPM1.404E-001 2.964E-001Max-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-001=-0.125 ，對(duì)應(yīng) c*x 以下的值，同時(shí)返回最 優(yōu)解：x1=1.9583,x2=2.0833, 發(fā)現(xiàn) x 確實(shí)有解，因?yàn)槠涿恳粋€(gè)元素都是非負(fù)的，而 且Ax=b,可以用命令min(x)和nor

5、m(Ax-b)來檢驗(yàn)。當(dāng)然，也會(huì)出現(xiàn)一些舍入誤 差，如下：norm(A*x-b)ans =1.7764e-015norm(A*(24*x)-24*b)ans = 0二次型和半定限制( Quadratic and semi definite constraint)在 sedumi 中,有可能強(qiáng)制加入二次型或者半定限制條件,即通過限制變量進(jìn) 入一個(gè)二次型和核或者半正定矩陣的核,這樣一個(gè)限制條件代替了線性規(guī)劃中的非 負(fù)性條件，需要x屬于K, 一個(gè)對(duì)稱核中，它是一個(gè)非負(fù)象限，二次型核和半正定 矩陣核的笛卡兒積。最優(yōu)化問題的標(biāo)準(zhǔn)形式為：minimizesuch that At-b（ （5） ）對(duì)偶形式為

6、：maximizeTFsuch （hat c A y XI.二次型核二次型核由以下形式來確定：QconefUhXiJe J? X 臚“訕 歸訕,考慮以下問題：min * Vi +鮒血I L P曲 啟Y1 *如卩，其中P為給定矩陣，q為給定向量，以上是魯棒最小均方問題。其中決策變 量為標(biāo)量y1和y2，還有向量y3.這個(gè)問題有兩個(gè)二次型限制：如、q - Pgm） ） E Qcone.（略 ）電 0門就給定P和q，以下MATLA函數(shù)（rls.m ）將該問題表述成標(biāo)準(zhǔn)對(duì)偶問題。 陣被表述為轉(zhuǎn)置方向，即表述為 At。Rls.m% At,b,c,KI = rls（P,q）% Creates dual st

7、a ndard form for robust least squares problem Pu=q.function 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);

8、 .zeros(1,2+n); . zeros(n,2),-eye(n);c = sparse(c;0;1;zeros(n,1);K.q= K.q, 2+n; 我們可以注意到以上的方程中運(yùn)用的是系數(shù)數(shù)據(jù)存儲(chǔ)形式，為的是節(jié)省內(nèi)存。此外，還定義了一個(gè)結(jié)構(gòu)體 K,其中K.q列出了二次型核的維數(shù)。K結(jié)構(gòu)體將 被運(yùn)用于告知 SeDuM：i c-ATy 的元素并不被限制為非負(fù)當(dāng)他們都被用于線性規(guī)劃it : b*y cg cgprec中。反而言之，第一個(gè)行 K.q(1) 被限制在一個(gè)二次型核中，最后一行 K.q(2) 被限 制在另外一個(gè)二次型核中，這就是我們?cè)谥敖?duì)稱核的方法，而且建立了兩個(gè) 二次型限制

9、。作為數(shù)值仿真的例子，我們來來解決一個(gè)魯棒最小平方值的問題，其依靠 P 中的列 P = 3 1 4;0 1 1;-2 5 3; 1 4 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 =

10、 5, dim = 11, blocks = 3nnz(A) = 16 + 0, nnz(ADA) = 23, nnz(L) = 14gap delta rate t/tP* t/tD* feas4.00E+000 0.0001: -3.10E+000 3.01E-001 0.000 0.0751 0.9900 0.9900 1.55 1 1 3.4E-0012 : -3.33E+000 6.02E-003 0.000 0.0200 0.99000.9900 1.03 1 1 6.2E-0033 : -3.33E+000 1.31E-005 0.000 0.0022 0.99900.9990

11、1.00 1 1 1.3E-0054 : -3.33E+000 1.26E-006 0.367 0.0958 0.99000.9900 1.00 1 1 1.3E-0065 : -3.33E+000 3.81E-008 0.000 0.0303 0.99000.9900 1.00 1 1 3.8E-0086 : -3.33E+000 2.24E-009 0.425 0.0587 0.99020.9900 1.00 1 1 2.3E-009iter seconds digits c*x b*y6 0.0 9.4 -3.3329085951e+000 -3.3329085963e+000|Ax-b

12、| =9.5e-010, Ay-c_+ =2.1E-009, |x|= 2.0e+000, |y|=2.5e+000Detailed timing (sec)Pre IPM Post3.120E-002 3.120E-002 0.000E+000Max-norms: |b|=1, |c| = 3,Cholesky |add|=0, |skip| = 0, |L.L| = 1.在以上這些SeDuMi的語句中，我們發(fā)現(xiàn)了一個(gè)新的輸入量 viz.K，這個(gè)量使 得SeDuMi能夠解決一個(gè)（5）和（6）形式的最優(yōu)化問題，其中對(duì)稱核 K被描述為 結(jié)構(gòu)體K。沒有K的話，SeDuMi只能解決（1）和（2）的問

13、題。我們可以通過驗(yàn)證（ 7）中的不等式來確定結(jié)果 y 是否能夠滿足條件（ 9）。 然而，SeDuM提供了一種更簡單的方法eigK，這個(gè)函數(shù)返回相對(duì)應(yīng)于對(duì)稱核的矩 陣的特征值。一個(gè)對(duì)稱形核包括哪些僅有非負(fù)特征值的向量，其中參見 Faraut and Koranyi9 。如下，我們能因此檢驗(yàn)可行性和最優(yōu)性：eigK（x,K）,eigK（c-At*y,K）ans =0.0000 -0.00001.41423.23070.0000 -0.00001.41421.4827x*（c-At*y）ans =3.4744e-009對(duì)于對(duì)稱核K來說，總是有xTz=0,對(duì)所有屬于K的z和x成立，因此，x 提供里一種

14、在線性規(guī)劃中對(duì)y最優(yōu)化的認(rèn)證。Farkas的對(duì)偶方案的分解也是采用 同樣的思想。具體參見15的方法。然而，一個(gè)奇怪的現(xiàn)象出現(xiàn)了，x和y幾乎是 可解得，然而cTx-bTy大多是負(fù)的(|x|和或者|y|然后肯定會(huì)是非常大的)， 根據(jù)原理，SeDuM將會(huì)返回一個(gè)無窮大的精確的 digits。這個(gè)現(xiàn)象在14和23 中也有解釋。我們需要讓這個(gè)最優(yōu)化模型有著非負(fù)的和二次型的核限制，而這一點(diǎn)是可行 的。比如，我們可以以上例子中的限制y31 ；血：巧+X is psd從幾何學(xué)上來說，Rcone僅僅是Qcone的轉(zhuǎn)置。Rcone的這種特殊形式對(duì)于建 立凸面體二次型方程十分方便。也就是說，將現(xiàn)行等式限制“x1=1

15、”加入到模型中，我們得到限制：通過該模型，我們能用x2作為的緊上界。分?jǐn)?shù)同樣也是由Rcone方便作為限制。舉例說明，我們可以最小化 1/x1對(duì)于x10來解決模型 J. Ji + J； 0.注意到該問題是無解的：1/x1下確界是0，對(duì)x1趨近與無窮大來說：clear K;c= 0,1,0；b = sqrt(2); A = 0, 0, 1; K.r = 3;x,y,i nfo=sedumi(A,b,c,K);x(2),x(1)*x (2)ans =6.5278e-005ans =1.0695你可能會(huì)發(fā)現(xiàn)x2根本沒有足夠趨近于0,而且x1也不是等于無窮大。然 而，這個(gè)原始解是可解的，對(duì)偶解幾乎是可解

16、的，而且對(duì)偶gap幾乎是負(fù)的。這就解釋了一個(gè)誤差范圍困難，而且對(duì)這種不規(guī)則的問題是很常見的。在Section 5中，我們可以看看到底如何獲得一個(gè)更精確的解決方案，通過設(shè)計(jì)一個(gè)選項(xiàng)參數(shù)，pars.eps上述例子可以解釋，K.r來列出Rcone限制的維數(shù)，類似于 Qcone限制中K.q 的定義，設(shè)定了 K.l，K.q，K.r也就引出了以下形式的對(duì)稱核。fC -胖x (Qcone x x Qcone) x (Rcone x - * x Reone).舉例說明，我們可以添加界 x1X=vec(1,5,-3;5,2,-9;-3,-9,4)X =15-352-9-3-94Vec的逆運(yùn)算就是mat。因此，如

17、果x是一個(gè)n2長度的向量，而mat (x) 就創(chuàng)建了一個(gè)nXn的矩陣，然后依次填入x向量的元素，比如：mat(X)ans =15-352-9-3-94以下MATLA函數(shù)產(chǎn)生了一個(gè)問題(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) -c = vec(spdiags(m-1:-l:0),0

18、,m,m);% - Let e be all-1, and allocate space for the A-matrix -e = ones(m,1);At = sparse( , ,mA2,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告訴SeDuM我們需要一個(gè)mXm的矩陣mat (x)，而且滿足對(duì)稱 正定的。我們能解決如下的問題：b = 2; 0.2; -0.3;At, b, c ,K = specfac(b) ;iter seconds

19、digitsc*xb*yfeasx,y,info = sedumi(At,b,c,K);得到的結(jié)果為：it : b*y gap delta rate t/tP* t/tD* cg cgprec0 : 1.69E+001 0.0001 : -2.44E-001 3.91E+000 0.000 0.2312 0.90000.9000 1.39 1 1 2.0E+0002 : 1.86E-001 3.22E-001 0.000 0.0824 0.99000.9900 1.31 1 1 4.0E-0013 : 1.23E-001 1.12E-004 0.000 0.0003 0.99990.9999

20、0.99 1 1 2.6E-0044 : 1.23E-001 5.18E-006 0.000 0.0463 0.99000.9900 1.00 1 1 1.2E-0055 : 1.23E-001 2.36E-007 0.000 0.0455 0.99000.9900 1.00 1 1 6.8E-0076 : 1.23E-001 1.90E-008 0.327 0.0806 0.99000.9900 1.00 2 2 5.5E-0087 : 1.23E-001 8.16E-010 0.000 0.0429 0.99060.9900 0.99 2 2 2.2E-0091.4E-010, |x|=

21、2.0e+000, |y|=Post3.120E-0020.2 8.5 1.2273256559e-001 1.2273256524e-001|Ax-b| = 1.1e-010, Ay-c_+ = 7.6e-001Detailed timing (sec)Pre IPM3.120E-002 2.028E-001Max-norms: |b|=2, |c| = 2,Cholesky |add|=0, |skip| = 0, |L.L| = 1.為檢驗(yàn)正定性，可以采用 MATLAB的eig或者SeDuMi中的eigK: eig(mat(x),eigK(x,K)ans =0.00000.00000.00000.00002.00002.0000其中 eigK 更加方便，特別是當(dāng)有多個(gè)正定條件限制的時(shí)候，或者也有非負(fù)和二次型核限制。SeDuM將總是產(chǎn)生對(duì)稱決策變量，比如 mat (x)就是對(duì)稱的。女口 果不明確地加入對(duì)稱性限制，比如 xij -xji=0 。最多，這些限制將會(huì)被 SeDuMi 在 A 矩陣中