這一段被老板逼着論文開題，自己找方向比較着急，最后選擇了供應鏈控制理論的一個方向。我要寫的論文，用到了Matlab的LMI工具，以及某篇論文中的H-inf穩定定理。自己好好研究了好長時間，怎么也無法實現該論文當中的算例。研究了一個多月，自己簡直就快崩潰了，也搞不定問題。我很是懷疑自己的選題是不是正確，並且懷疑自己是不是選的太難了。如果連論文中的算例都無法實現，如何把該模型應用到自己論文當中呢？功夫不負有心人，昨日我加入了Mathworks的Matlab的Newsgroup，結果碰見一牛人Johan，立即就把論文中的算例給寫成程序。但是他做出的結果和論文仍然有差別，我仍有點不甘心，人家的論文發表在Automatica上，如果連這種期刊都水的要命，那么就沒有什么學術水平了。
今天中午，仍然不甘心，老爸給我打了電話讓我看紅場閱兵，於是我邊看PPMate邊漫無邊際的核對着自己的程序。終於做出了和算例一致的結果。
我搜出來的都是一些簡單的算例，並且機會沒有中文教程，我在這里就斗膽把自己的體會寫出來，試着給大家提供一點參考。
LMI：Linear Matrix Inequality，就是線性矩陣不等式。
在Matlab當中，我們可以采用圖形界面的lmiedit命令，來調用GUI接口，但是我認為采用程序的方式更方便（也因為我不懂這個lmiedit的GUI）。
對於LMI Lab， 其中有三種求解器（solver）： feasp，mincx和gevp。
每個求解器針對不同的問題：
feasp：解決可行性問題（feasibility problem），例如：A(x)<b(x)。< font="">
mincx：在線性矩陣不等式的限制下解決最小化問題（Minimization of a linear objective under LMI constraints），例如最小化c'x，在限制條件A(x) < B(x)下。
gevp：解決廣義特征值最小化問題。例如：最小化lambda，在0<b(x),a(x)<lamba*b(x)限制條件下。< font="">
要解決一個LMI問題，首要的就是要把線性矩陣不等式表示出來。
對於以下類型的任意的LMI問題
N' * L(X1, . . . , XK) * N < M' * R(X1, . . . , XK) * M
其中X1, . . . , XK是結構已經事先確定的矩陣變量。左側和右側的外部因子（outer factors）N和M是給定的具有相同維數的矩陣。
左側和右側的內部因子（inner factors）L(.)和R(.)是具有相同結構的對稱塊矩陣。每一個塊由X1, . . . , XK以及它們的轉置組合而成形成的。
解決LMI問題的步驟有兩個：
1、定義維數以及每一個矩陣的結構，也就是定義X1, . . . , XK。
2、描述每一個LMI的每一項內容（Describe the term content of each LMI）
此處介紹兩個術語：
矩陣變量（Matrix Variables）：例如你要求解X滿足A(x)<b(x)，那么x就叫做矩陣< font="">
變量。
項（Terms）：項是常量或者變量（Terms are either constant or variable）。
常項（Constant Terms）是確定的矩陣。可變項（Variable Terms)是哪些含有矩陣變
量的項，例如：X*A, X*C'。如果是X*A + X*C'，那么記得要把它當成兩項來處理。
好了廢話不說了，讓我們來看個例子吧（下面是一線性時滯系統）。
語法：X = lmivar(type,struct)
type=1: 定義塊對角的對稱矩陣。 每一個對角塊或者是全矩陣<任意對稱矩陣>，標量<單位矩陣的乘積>，或者是零陣。
如果X有R個對角塊，那么后面這個struct就應該是一個Rx2階的的矩陣，在此矩陣中，struct(r,1)表示第r個塊的大小，struct(r,2) 表示第r個塊的類型<1--全矩陣，0--標量，-1--零陣）。
比如一個矩陣有兩個對角塊，其中一個是2x2的全對稱矩陣，第二個是1x1的一個標量，那么該矩陣變量應該表示為X = lmivar(1, [2 1; 1 0]) 。
type=2: mxn階的矩陣，只需要寫作struct = [m,n]即可。
type=3: 其它類型。針對類型3，X的每一個條目（each entry of X）被定義為0或者是+(-)xn，此處xn代表了第n個決策變量。
那么針對我們的例子，我們如此定義變量：
% Q is a symmetric matrix, has a block size of 2 and this block is symmetric
Q = lmivar(1, [2 1]);
% S1 a symmeric matrix, size 2
S1 = lmivar(1, [2 1]);
% S2 is 1 by 1 matrix
S2 = lmivar(1, [1 0]);
% Type of 2, size 1 by 2
M = lmivar(2, [1 2]);
定義完成變量之后，我們就該用lmiterm來描述LMI中的每一個項了。Matlab的官方文檔提示我們，如果要描述一個LMI只需要描述上三角或者下三角元素就可以了，否則會描述成另一個LMI。
When describing an LMI with several blocks, remember to specify only the terms in the blocks on or below the diagonal (or equivalently, only the terms in blocks on or above the diagonal).
語法為：lmiterm(termID,A,B,flag)
termID是一個四維整數向量，來表示該項的位置和包含了哪些矩陣變量。
termID(1)可以為+p或者-p，+p代表了這個項位於第p個線性矩陣不等式的左邊，-p代表了這個項位於第p個線性矩陣不等式的右邊。注意：按照慣例來講，左邊通常指較小的那邊。
termID(2:3):
1、對於外部變量來說，取值為[0,0];
2、對於左邊或者右邊的內部變量來說，如果該項在(i,j)位置，取值[i,j]
termID(4):
1、對於外部變量，取值為0
2、對於A*X*B，取值X
3、對於A*X'*B，取值-X
flag(可選，值為s）:
因為：(A*X*B) + (A*X*B)T = A*X*B + B'*X'*A'，所以采用s來進行簡寫。
比如：針對A*X + X'*A'
我們采用笨方法：
lmiterm([1 1 1 X],A,1)
lmiterm([1 1 1 -X],1,A')
那么簡寫就是lmiterm([1 1 1 X],A,1,'s')
接下來我們就看該論文中的算例吧：(1,1)位置是
-Q+Bd*S2*Bd'+Ad*S1*Ad';
我們應該表示為：
% pos in (1, 1)
lmiterm([1 1 1 Q], -1, 1);
lmiterm([1 1 1 S2], Bd, Bd');
lmiterm([1 1 1 S1], Ad, Ad');
其它位置仿照寫就行了，不懂了多看幫助文檔。
把每一個項都定義以后，要記得
lmis = getlmis;
[tmin, feas] = feasp(lmis)
getlmis:是在完成定義變量和項之后，LMI系統的內部表示就可以通過此命令獲得（After completing the description of a given LMI system with lmivar and lmiterm, its internal representation lmisys is obtained with the command）。
feasp是調用feasp求解器，看有沒有可行解。feas就是可行解。
下面我把代碼貼上去，那些常數矩陣都在此源程序中定義了。
A = [2 1; 0 1];
Ad = [0.2 0.1; 0 0.1];
B1 = [0.1 0.1]';
B2 = [1 1]';
Bd = [0.1 0.1]';
C = [1, 1];
Cd = [0.1, 0.1];
D11 = 0.1;
D12 = 1;
Dd = 0.1;
gammar = 1;
% Initial a LMI system
setlmis([]);
% Define Variables
% Q is a symmetric matrix, has a block size of 2 and this block is symmetric
Q = lmivar(1, [2 1]);
% S1 a symmeric matrix, size 2
S1 = lmivar(1, [2 1]);
% S2 is 1 by 1 matrix
S2 = lmivar(1, [1 0]);
% Type of 2, size 1 by 2
M = lmivar(2, [1 2]);
% Q, S1, S2 > 0
lmiterm([-2 1 1 Q], 1, 1);
lmiterm([-3 1 1 S1], 1, 1);
lmiterm([-4 1 1 S2], 1, 1);
% pos in (1, 1)
lmiterm([1 1 1 Q], -1, 1);
lmiterm([1 1 1 S2], Bd, Bd');
lmiterm([1 1 1 S1], Ad, Ad');
% pos (1, 2)
lmiterm([1 1 2 Q], A, 1);
lmiterm([1 1 2 M], B2, 1);
% pos(1, 3)
lmiterm([1 1 3 0], B1);
% pos(1, 4)
lmiterm([1 1 4 S2], Bd, Dd');
lmiterm([1 1 4 S1], Ad, Cd');
% pos(2, 2)
lmiterm([1 2 2 Q], -1, 1);
% pos(2, 4)
lmiterm([1 2 4 Q], 1, C');
lmiterm([1 2 4 -M], 1, D12');
% pos(2, 5)
lmiterm([1 2 5 -M], 1, 1);
% pos(2, 6)
lmiterm([1 2 6 Q], 1, 1);
% pos(3, 3)
lmiterm([1 3 3 0], -(gammar^2));
% pos(3, 4)
lmiterm([1 3 4 0], D11');
% pos(4, 4)
lmiterm([1 4 4 0], -1);
lmiterm([1 4 4 S1], Cd, Cd');
lmiterm([1 4 4 S2], Dd, Dd');
lmiterm([1 5 5 S2], -1, 1);
lmiterm([1 6 6 S1], -1, 1);
lmis = getlmis;
[tmin, feas] = feasp(lmis)
運行后，就調用dec2mat把決策變量轉化為矩陣形式。
Q = dec2mat(lmis, feas, Q)
Q =
1.9253 -2.2338
-2.2338 9.1054
可以看到，和論文中的一樣。