SeDuMi详细教程_sedumi详细教程
SeDuMi详细教程由刀豆文库小编整理，希望给你工作、学习、生活带来方便，猜你可能喜欢“sedumi详细教程”。
SeDuMi User Guide
线性规划（LP问题：Linear Programming）
将需要解决的问题表述成原始问题：
或者对偶问题：
例子：
为了解决该LP问题，我们必须添加一些松弛变量，比如x3和x4，我们即可在MATLAB 中输入b和c向量以及A矩阵:
现在我们就可以通过sedumi指令来求解这个问题 >> sedumi(A,b,c)eqs m = 2, order n = 5, dim = 5, blocks = 1 nnz(A)= 6 + 0, nnz(ADA)= 4, nnz(L)= 3 it :
b*y
gap
delta rate
t/tP* t/tD*
feas cg cg prec
0 : 6.02E+000 0.000:-1.20E+000 1.84E+000 0.000 0.3055 0.9000 0.9000
1.86 1 1 2.1E+000:-6.94E-002 5.15E-001 0.000 0.2803 0.9000 0.9000
2.14 1 1 8.6E-001:-1.21E-001 2.72E-002 0.000 0.0528 0.9900 0.9900
1.16 1 1 4.2E-002:-1.25E-001 8.69E-006 0.000 0.0003 0.9999 0.9999
1.02 1 1
iter seconds digits
c*x
b*y
0.3
Inf
-1.2500000000e-001-1.2500000000e-001 |Ax-b| = 0.0e+000, [Ay-c]_+ = 1.7E-016, |x|= 2.9e+000, |y|= 2.8e-001
Detailed timing(sec)
Pre
IPM
Post 1.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 表明最优值为-1.2500000000e-001=-0.125，对应c*x以下的值，同时返回最优解：
x1=1.9583,x2=2.0833,发现x确实有解，因为其每一个元素都是非负的，而且Ax=b，可以用命令min(x)和norm(Ax-b)来检验。当然，也会出现一些舍入误差，如下：
>> norm(A*x-b)ans =1.7764e-015 >> norm(A*(24*x)-24*b)ans = 0 二次型和半定限制（Quadratic and semi definite constraint）
在sedumi中，有可能强制加入二次型或者半定限制条件，即通过限制变量进入一个二次型和核或者半正定矩阵的核，这样一个限制条件代替了线性规划中的非负性条件，需要x属于K，一个对称核中，它是一个非负象限,二次型核和半正定矩阵核的笛卡儿积。最优化问题的标准形式为：
对偶形式为：
二次型核
二次型核由以下形式来确定：
考虑以下问题：
其中P为给定矩阵，q为给定向量，以上是鲁棒最小均方问题。其中决策变量为标量y1和y2，还有向量y3.这个问题有两个二次型限制：
给定P和q，以下MATLAB函数（rls.m）将该问题表述成标准对偶问题。A矩阵被表述为转置方向，即表述为At。Rls.m % [At,b,c,KI = rls(P,q)
% Creates dual standard form for robust least squares problem “Pu=q”.function [At ,b,c,K]= rls(P,q)[m, n] = size(P);
%--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)];c = sparse([c;0;1;zeros(n,1)]);K.q= [K.q, 2+n];
我们可以注意到以上的方程中运用的是系数数据存储形式，为的是节省内存。此外，还定义了一个结构体K，其中K.q列出了二次型核的维数。K结构体将被运用于告知SeDuMi：c-ATy的元素并不被限制为非负当他们都被用于线性规划中。反而言之，第一个行K.q(1)被限制在一个二次型核中，最后一行K.q(2)被限制在另外一个二次型核中，这就是我们在之前建立对称核的方法，而且建立了两个二次型限制。
作为数值仿真的例子，我们来来解决一个鲁棒最小平方值的问题，其依靠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);运行后出现的结果是
SeDuMi 1.1R3 by AdvOL, 2006 and Jos F.Sturm, 1998-2003.Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500 eqs m = 5, order n = 5, dim = 11, blocks = 3 nnz(A)= 16 + 0, nnz(ADA)= 23, nnz(L)= 14 it :
b*y
gap
delta rate
t/tP* t/tD*
feas cg cg prec
0 :
4.00E+000 0.000:-3.10E+000 3.01E-001 0.000 0.0751 0.9900 0.9900
1.55 1 1 3.4E-001:-3.33E+000 6.02E-003 0.000 0.0200 0.9900 0.9900
1.03 1 1 6.2E-003:-3.33E+000 1.31E-005 0.000 0.0022 0.9990 0.9990
1.00 1 1 1.3E-005:-3.33E+000 1.26E-006 0.367 0.0958 0.9900 0.9900
1.00 1 1 1.3E-006:-3.33E+000 3.81E-008 0.000 0.0303 0.9900 0.9900
1.00 1 1 3.8E-008:-3.33E+000 2.24E-009 0.425 0.0587 0.9902 0.9900
1.00 1 1 2.3E-009 iter seconds digits
c*x
b*y
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|= 2.5e+000
Detailed timing(sec)
Pre
IPM
Post 3.120E-002
3.120E-002
0.000E+000
Max-norms: ||b||=1, ||c|| = 3, Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.在以上这些SeDuMi的语句中，我们发现了一个新的输入量viz.K，这个量使得SeDuMi能够解决一个（5）和（6）形式的最优化问题，其中对称核K被描述为结构体K。没有K的话，SeDuMi只能解决（1）和（2）的问题。
我们可以通过验证（7）中的不等式来确定结果y是否能够满足条件（9）。然而，SeDuMi提供了一种更简单的方法eigK，这个函数返回相对应于对称核的矩阵的特征值。一个对称形核包括哪些仅有非负特征值的向量，其中参见Faraut and Koranyi[9]。
如下，我们能因此检验可行性和最优性： >> [ eigK(x, K), eigK(c-At*y, K)] ans =
0.0000
-0.0000
1.4142
3.2307
0.0000
-0.0000
1.4142
1.4827 >> x'*(c-At*y)ans =
3.4744e-009 对于对称核K来说，总是有xTz>=0，对所有属于K的z和x成立，因此，x提供里一种在线性规划中对y最优化的认证。Farkas的对偶方案的分解也是采用同样的思想。具体参见[15]的方法。然而，一个奇怪的现象出现了，x和y几乎是可解得，然而cTx-bTy大多是负的（||x||和或者||y||然后肯定会是非常大的），根据原理，SeDuMi将会返回一个无穷大的精确的digits。这个现象在[14]和[23]中也有解释。
我们需要让这个最优化模型有着非负的和二次型的核限制，而这一点是可行的。比如，我们可以以上例子中的限制y3[1]
a1 = zeros(1,length(y));a1(3)= 1;c= [-0.1;c];At = [a1;At];K.l=1;
[x,y,info] = sedumi(At,b,c,K);eigK(c-At*y,K)'
ans =
0.0000
0.0000
3.2307
-0.0000
1.4904 其中K.l是非负变量的数目，在这种情况下是1，按照惯例，非负的变量一般是第一个元素，因此我们的例子中K=R+*Qcone*Qcone，从eigK的输出中可以见到，该核有5个特征值：1是非负性限制，2是对每个二次型限制。我们需要说明K是一个5阶的对称核。（SeDuMi返回‘order= 6’，因为它们的自身的自对偶变化）
SeDuMi提供了另外一种二次型核的形式：
方便。也就是说，将现行等式限制“x1=1”加入到模型中，我们得到限制：
从几何学上来说，Rcone仅仅是Qcone的转置。Rcone的这种特殊形式对于建立凸面体二次型方程十分
通过该模型，我们能用x2作为的紧上界。分数同样也是由Rcone方便作为限制。举例说明，我们可以最小化1/x1对于x1>0来解决模型
注意到该问题是无解的：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 ans = 1.0695 你可能会发现x2根本没有足够趋近于0，而且x1也不是等于无穷大。然而，这个原始解是可解的，对偶解几乎是可解的，而且对偶gap几乎是负的。这就解释了一个误差范围困难，而且对这种不规则的问题是很常见的。在Section 5 中，我们可以看看到底如何获得一个更精确的解决方案，通过设计一个选项参数，pars.eps。
上述例子可以解释，K.r来列出Rcone限制的维数，类似于Qcone限制中K.q的定义，设定了K.l，K.q，K.r也就引出了以下形式的对称核。
举例说明，我们可以添加界‘x1
[x,y,info]=sedumi(A,b,c,K);
半正定核（the positive semidefinite cone）
半正定限制条件是SeDuMi建模的一个重要的限制集合，举个例子，考虑到以下问题：
其中，X是m×m对称矩阵，xij代表了第i行第j列的元素。长度m的向量b被给定。缩写‘psd’表示‘positive semidefinite’，以上最优化问题变为一个自相关向量b的最小相位的谱分解问题，参见[4]。鉴于SeDuMi善于处理决策变量向量的问题，问题（11）以一个m×m对称矩阵的决策变量的方式被提出。这个小问题可以用有名的向量化方法来解决。向量化方法可以通过vec和mat函数实现，这也是SeDuMi的一部分。函数vec（X）创建了一个长向量，通过堆叠了矩阵X的每一列，比如： >> X=vec([1,5,-3;5,2,-9;-3,-9,4])' X =
-9Vec的逆运算就是mat。因此，如果x是一个n^2长度的向量，而mat（x）就创建了一个n×n的矩阵，然后依次填入x向量的元素，比如： >> mat(X)ans =
-9-3
-9以下MATLAB函数产生了一个问题（11）的标准的原始形式： % [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,m,m));
%-----Let e be all-1, and allocate space for the A-matrix-------e = ones(m,1);
At = sparse([], [], [] ,m^2,m,m*(m+l)/2);%-------sum(diag(): k))= b(k)-------------for k= 1:m
At(:,k)= vec(spdiags(e,k-l,m,m));end K.s = m;字段K.s=m告诉SeDuMi我们需要一个m×m的矩阵mat（x），而且满足对称正定的。我们能解决如下的问题：
b = [2;0.2;-0.3];[At, b, c ,K] = specfac(b);[x,y,info] = sedumi(At,b,c,K);得到的结果为：
it :
b*y
gap
delta rate
t/tP* t/tD*
feas cg cg prec
0 :
1.69E+001 0.000:-2.44E-001 3.91E+000 0.000 0.2312 0.9000 0.9000
1.39 1 1 2.0E+000: 1.86E-001 3.22E-001 0.000 0.0824 0.9900 0.9900
1.31 1 1 4.0E-001: 1.23E-001 1.12E-004 0.000 0.0003 0.9999 0.9999
0.99 1 1 2.6E-004: 1.23E-001 5.18E-006 0.000 0.0463 0.9900 0.9900
1.00 1 1 1.2E-005: 1.23E-001 2.36E-007 0.000 0.0455 0.9900 0.9900
1.00 1 1 6.8E-007: 1.23E-001 1.90E-008 0.327 0.0806 0.9900 0.9900
1.00 2 2 5.5E-008: 1.23E-001 8.16E-010 0.000 0.0429 0.9906 0.9900
0.99 2 2 2.2E-009
iter seconds digits
c*x
b*y
0.2
8.5 1.2273256559e-001 1.2273256524e-001 |Ax-b| = 1.1e-010, [Ay-c]_+ = 1.4E-010, |x|= 2.0e+000, |y|= 7.6e-001
Detailed timing(sec)
Pre
IPM
Post 3.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.0000
0.0000
0.0000 2.0000
2.0000 其中eigK更加方便，特别是当有多个正定条件限制的时候，或者也有非负和二次型核限制。SeDuMI将总是产生对称决策变量，比如mat（x）就是对称的。如果不明确地加入对称性限制，比如‘xij-xji=0’。最多，这些限制将会被SeDuMi在A矩阵中剔除。
然而，对偶问题的解c-ATy不需要对称，就像我们处理的数字例子中： >> mat(c-At*y)ans =
2.0727
-0.3130
0.6849
0
1.0727
-0.3130
0
0
0.0727 在这种情况下，对偶问题是一个上三角的形式，因为mat（c）是对角线的，而且mat（At（：，k））是对于k=1,2…m是上三角的。我们让Z=mat（c-ATy），SeDuMi限制了Z的对称部分，是（Z+ZT）/2成为正定的。函数eigK变成这个对称部分的特征根，因此, >> eigK(c-At*y,K)ans =
-0.0000
1.0583 2.1597 产生与以下一样的结果: >> Z=mat(c-At*y);eig((Z+Z')/2)ans =
-0.0000
1.0583 2.1597 我们可以注意到这个问题：
也就是说，for k=1:size(At,2),Ak=mat(At(:,k));
At(:,k)=vec(Ak+Ak')/2;
end 然后可以通过SeDuMi产生一样的解x和y。然而因为A矩阵限制是对称的，我们现在发现mat（c-At*y）也是对称的，即为(Z+Z')/2。
对于有多个正定限制的情况，K.s可以列出相关矩阵的阶数，类似于K.q和K.r的定义。正定决策变量常会是x和c-ATy的最后的元素。，因为X是对称的。因此，我们可以如下改变A矩阵：
其中Scone表示正定矩阵的核。
