LMI工具箱介绍
线性矩阵不等式(LMI)工具箱是求解一般线性矩阵不等式问题的一个高性能软件
包。由于其面向结构的线性矩阵不等式
示方式,使得各种线性矩阵不等式能够以自然块
矩阵的形式加以描述。一个线性矩阵不等式问题一旦确定,就可以通过调用适当的线性矩
阵不等式求解器来对这个问题进行数值求解。
LMI 工具箱提供了确定、处理和数值求解线性矩阵不等式的一些工具,它们主要用
于:
z 以自然块矩阵形式来直接描述线性矩阵不等式;
z 获取关于现有的线性矩阵不等式系统的信息;
z 修改现有的线性矩阵不等式系统;
z 求解三个一般的线性矩阵不等式问题;
z 验证结果。
本附录将详细介绍 LMI 工具箱提供的用于解决以上各个问题的相关函数和命令。
A.1 线性矩阵不等式及相关术语
一个线性矩阵不等式就是具有以下一般形式的一个矩阵不等式:
0<+++= NNxx LLLxL "110)( (1)
其中: 是给定的对称常数矩阵, 是未知变量,称为决策变量,NLLL ,,, 10 " Nxx ,,1 "
∈= T1 ],,[ Nxx "x NR 是由决策变量构成的向量,称为决策向量。
尽管表达式(1)是线性矩阵不等式的一个一般表示式,但在大多数实际应用中,线
性矩阵不等式常常不是以一般表示式(1)的形式出现,而是具有以下形式:
),,(),,( 11 nn XXRXXL "" <
其中的 和 是矩阵变量 的仿射函数,通过适当的代数运算,上式可以写
成线性矩阵不等式的一般表示式(1)的形式。例如,在系统稳定性问题中经常遇到的
Lyapunov 矩阵不等式
)(⋅L )(⋅R nXX ,,1 "
0<+ XAXAT (2)
也是一个线性矩阵不等式,其中的 是一个矩阵变量。我们以一个二阶矩阵
为例,将矩阵不等式(2)写成一般表示式(1)的形式。针对二阶矩阵不
X
⎥⎦
⎤⎢⎣
⎡
−
−=
20
21
A
等式(2),对应的矩阵变量 是一个二阶的对称矩阵, ,不等式(2)中
的决策变量是矩阵 中的独立元 。根据对策性,矩阵变量 可以写成
X ⎥⎦
⎤⎢⎣
⎡=
32
21
xx
xx
X
X 321 xxx 、、 X
⎥⎦
⎤⎢⎣
⎡+⎥⎦
⎤⎢⎣
⎡+⎥⎦
⎤⎢⎣
⎡=
10
00
01
10
00
01
321 xxxX
将矩阵 A和上式代入矩阵不等式(2),经整理,可得
0<⎥⎦
⎤⎢⎣
⎡
−+⎥⎦
⎤⎢⎣
⎡
−
−+⎥⎦
⎤⎢⎣
⎡−
40
00
43
30
02
22
321 xxx (3)
这样就将矩阵不等式(2)写成了线性矩阵不等式的表示式(1)。显然,与 Lyapunov 矩
阵不等式(2)相比,表示式(3)缺少了许多控制中的直观意义。另外,(3)式涉及到
的矩阵也比(2)式中的多。如果矩阵 是 阶的,则(3)式中的系数矩阵一般有A n
2)1( +nn 个。因此,这样的表达式在计算机中将占用更多的存储空间。由于这样的一些
原因,LMI 工具箱中的函数采用线性矩阵不等式的结构表示。例如,Lyapunov 矩阵不等
式(2)就以矩阵变量 的不等式来表示,而不是用其一般形式(3)来表示。 X
一般的,一个线性矩阵不等式具有块矩阵的形式,其中每一个块都是矩阵变量的仿射
函数。以下通过一个例子来说明有关描述一个线性矩阵不等式的术语。
考虑 控制中的一个线性矩阵不等式: ∞H
0<
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−
−
+
N
IDB
DICX
BXCXAXA
N
γ
γ
TT
TT
T
其中: 是给定的矩阵,NDCBA 、、、、 ∈= TXX nn×R 和 R∈γ 是问题的变量。
z N称为外因子,块矩阵
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−
−
+
=
IDB
DICX
BXCXAXA
XL
γ
γγ
TT
TT
),(
称为内因子。外因子可以不是一个正方矩阵,它在许多问题中常常不出现。
z X 和γ 是问题的矩阵变量。注意标量也可以看成是一个 11× 维的矩阵。
z 内因子 ),( γXL 是一个对称块矩阵。根据对称性, ),( γXL 可以由对角线及其上方
的块矩阵完全确定。
z ),( γXL 中的每一块都是矩阵变量 X 和 γ 的仿射函数。这一函数由常数项和变量
项这两类基本项组成,其中常数项就是常数矩阵或以一些常数矩阵组成的算术表
达式,例如 ),( γXL 中的 和B D ;变量项是包含一个矩阵变量的项,例如
IXA γ−, 等。
一个线性矩阵不等式不论多么复杂,都可以通过描述其中每一块的各项内容来确定这
个线性矩阵不等式。
A.2 线性矩阵不等式的确定
LMI 工具箱可以处理具有以下一般形式的线性矩阵不等式:
MXXRMNXXLN ),,(),,( 1
T
1
T
KK "" <
其中: 是具有一定结构的矩阵变量,左、右外因子KXX ,,1 " N和M 是具有相同维数的
给定矩阵,左、右内因子 和)(⋅L )(⋅R 是具有相同块结构的对称块矩阵。
注意在线性矩阵不等式的描述中,左边总是指不等式较小的一边,例如对线性矩阵不
等式 ,0>X X 称为是不等式的右边,0 称为是不等式的左边,常表示成 。 X<0
要确定一个线性矩阵不等式系统,需要做以下两步:
1.给出每个矩阵变量 的维数和结构; KXX ,,1 "
2.描述每一个线性矩阵不等式中各个项的内容。
这个过程产生所描述线性矩阵不等式系统的一个内部表示,它以一个单一向量的形式
储存在计算机内,通常用一个名字,例如 lmisys 来表示。该内部表示 lmisys 可以在后面处
理这个线性矩阵不等式时调用。
下面将通过 LMI 工具箱中的一个例子来说明线性矩阵不等式系统的确定。运行
lmidem可以看到这个例子的完整描述。
例 1:考虑一个具有 4 个输入、4 个输出和 6 个状态的稳定传递函数
BAICG 1)()( −−= ss (4)
和一组具有以下块对角结构的输入/输出尺度矩阵 D:
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
=
54
32
1
1
00
00
000
000
dd
dd
d
d
D (5)
则在具有时变不确定性系统的鲁棒稳定性
中提出了以下问题:
寻找一个具有结构(5)的尺度矩阵 D,使得 1)(sup 1 <−DDG ω
ω
j 。
可以证明:这样一个问题可以转化成一个线性矩阵不等式系统的可行性问题,即寻找
两个对称矩阵 X 66R ×∈ 和 DDS T= 44×∈R ,使得
0<⎥⎦
⎤⎢⎣
⎡
−
++
SXB
XBSCCXAXA
T
TT
(6)
0>X (7)
IS > (8)
用命令 lmivar 和 lmiterm 给出线性矩阵不等式系统(6)~(8)的内部描述如下:
setlmis([])
X=lmivar(1,[6 1])
S=lmivar(1,[2 0;2 1])
% 1st LMI
lmiterm([1 1 1 X],1,A,’s’)
lmiterm([1 1 1 S],C’,C)
lmiterm([1 1 2 X],1,B)
lmiterm([1 2 2 S],-1,1)
% 2nd LMI
lmiterm([-2 1 1 X],1,1)
% 3rd LMI
lmiterm([-3 1 1 S],1,1)
lmiterm([3 1 1 0],1)
lmisys=getlmis
其中:函数 lmivar 定义了两个矩阵变量 X 和 ,lmiterm 则描述了每一个线性矩阵不等式
中各项的内容。getlmis 回到了这个线性矩阵不等式系统的内部表示 lmisys,lmisys 也称为
是储存在机器内部的线性矩阵不等式系统的名称。以下将详细介绍这几个函数的功能和用
法。
S
setlmis 和 getlmis
一个线性矩阵不等式系统的描述以 setlmis 开始,以 getlmis 结束。当要确定一个新的
系统时,输入:
setlmis([])
如果需要将一个线性矩阵不等式添加到一个名为 lmiso 的现有的线性矩阵不等式系统中,
则输入:
setlmis(lmiso)
当线性矩阵不等式系统被完全确定好后,输入:
lmisys=getlmis
该命令返回这个线性矩阵不等式系统的内部表示 lmisys。
lmivar
函数 lmivar 用来描述出现在线性矩阵不等式系统中的矩阵变量,每一次只能描述一个
矩阵变量。矩阵变量的描述包括该矩阵变量的结构。该函数的一般表达式是:
X=lmivar(type,struct)
这一函数定义了一个新的矩阵变量 X ,X 是该矩阵变量的变量名。函数中的第一个输
入量 type 确定了矩阵变量 X的类型,第二个输入量 struct 进一步根据变量 X 的类型给出该
变量的结构。变量的类型分成三类:
Type =1:对称块对角结构。这种结构对应于具有以下形式的矩阵变量:
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
rD
D
D
"
#%##
"
"
00
00
00
2
1
其中对角线上的每一个矩阵块 是方阵,它可以是零矩阵、对称矩阵或数量矩阵。这种
结构也包含了通常意义的对称矩阵和数量矩阵(分别相当于只有一块)。此时,struct 是
一个
jD
2×r 维的矩阵。如果该矩阵的第 i行是 ,则其中的 表示对称矩阵块 的阶
数,而 n只能取 1、0 或 ,其中
),( nm m iD
1− 1=n 表示 是一个满的对称矩阵(或无结构的对称矩
阵), 表示 是一个数量矩阵,
iD
0=n iD 1−=n 表示 是一个零矩阵。 iD
Type =2:长方型结构。这种结构对应于任意的长方矩阵。此时,struct= 表示矩
阵的维数。
),( nm
Type =3:其他结构。这种结构用来描述更加复杂的矩阵,也可以用于描述矩阵变量
之间的一些关联。 X 的每一个元或者是 0,或者是 nx± ,其中 是第 个决策变量。相应
的,struct 是一个和变量
nx n
X 有相同维数的矩阵,其中的每一个元取值如下:
struct
⎪⎩
⎪⎨
⎧
−=−
=
=
=
n
n
xjin
xjin
ji
ji
),(,
),(,
0),(,0
),(
X
X
X
如果
如果
如果
例 2:考虑具有三个矩阵变量 和 的线性矩阵不等式系统,其中 21 XX 、 3X
z 是一个1X 33× 维的对称矩阵;
z 是一个 维的长方矩阵; 2X 42×
z ,其中
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡Δ
=
22
13
I
X
δ
δ
00
00
00
Δ是 55× 维的对称矩阵, 1δ 和 2δ 是两个标量, 表示
维的单位矩阵。
2I
22×
可以应用 lmivar 来定义这些矩阵变量:
setlmis([])
X1=lmivar(1,[3 1])
X2=lmivar(2,[2 4])
X3=lmivar(1,[5 1;1 0;2 0])
lmiterm
在确定了矩阵变量之后,还需要确定每一个线性矩阵不等式中各项的内容。线性矩阵
不等式的项指构成这个线性矩阵不等式的块矩阵中的加项。这些项可以分成三类:
1.常数项;
2.变量项,即包含了矩阵变量的项,例如(3)式中的 XAT 和 。一般的变量
项具有形式 ,其中的
SCC T
PXQ X 是一个变量, P 和Q是给定的矩阵,分别称为该变量项的左
系数和右系数;
3.外因子。
在描述一个具有多个块的线性矩阵不等式时,LMI 工具箱提供了这样的功能,即只需
要确定对角线上和对角线上方的项的内容,或者只描述对角线上和对角线下方的项的内
容,其他部分项的内容可以根据线性矩阵不等式的对称性得到。
用命令 lmiterm 每次可以确定线性矩阵不等式的一个项的内容。例如,对线性矩阵不
等式
0<⎥⎦
⎤⎢⎣
⎡
−
++
SXB
XBSCCXAXA
T
TT
可以用以下一组命令来描述:
lmiterm([1 1 1 X],1,A,’s’)
lmiterm([1 1 1 S],C’,C)
lmiterm([1 1 2 X],1,B)
lmiterm([1 2 2 S],-1,1)
这些命令依次描述了项 和XBSCCXAXA 、、 TT + S− 。在每一条命令中,第 1 项是
一个四元向量,它刻画了所描述的项所在的位置和特征:
z 第 1 个元表示所描述的项属于哪一个线性矩阵不等式。值 m 表示第 m 个不等式的
左边,-m表示第 m个不等式的右边。
z 第 2 和第 3 个元表示所描述的项所在块的位置。例如,向量[1 1 2 1]表示所描述的
项位于第一个线性矩阵不等式左边内因子的块(1,2)中。第 2 和第 3 个元均取
零表示所描述的项在外因子中。
z 最后一个元表明了所描述的项是常数项还是变量项。如果是变量项,则进一步说
明涉及哪一个变量。0 表示常数项,k 表示所描述的项包含第 k 个矩阵变量 ,-
k 则表示包含矩阵变量 的转置 (在例 1 中,
kX
kX
T
kX X 是第 1 个变量, 是第 2 个
变量,它们按确定的先后顺序排列)。
S
lmiterm 的第 2 项和第 3 项包含了数据(常数项的值,外因子,变量项 或
中的左、右系数)。第 4 项是可选择的,且只能是's'。
PXQ QPX T
在描述项的内容时,有一些简化的方法。
1.零块可以省略描述;
2.可以通过在命令 lmiterm 中外加一个分量's',使得可以只用一个命令 lmiterm 就能
Administrator
高亮
描述一个变量项与该变量项的转置的和。例如,上面的第一个命令描述了 XAXA +T 。
3.可以用一个标量值来表示一个数量矩阵,即用α 表示数量矩阵 Iα ,其中α 是一个
标量。如例 1 中的第 3 个不等式 被描述成 IS >
lmiterm([-3 1 1 S],1,1)
lmiterm([3 1 1 0],1)
为了便于阅读,也可以用线性矩阵不等式和矩阵变量的名称来表示对应的线性矩阵不
等式和矩阵变量。矩阵变量的变量名可以用命令 lmivar 来赋值,线性矩阵不等式的名称则
可以用函数 newlmi 来确定。这些标识符可以用在命令 lmiterm 中以表示相应的线性矩阵不
等式或矩阵变量。对例 1 中的线性矩阵不等式系统,采用名称的相应描述如下:
setlmis([])
X=lmivar(1,[6 1])
S=lmivar(1,[2 0;2 1])
BRL=newlmi
lmiterm([BRL 1 1 X],1,A,’s’)
lmiterm([BRL 1 1 S],C’,C)
lmiterm([BRL 1 2 X],1,B)
lmiterm([BRL 2 2 S],-1,1)
Xpos=newlmi
lmiterm([-Xpos 1 1 X],1,1)
Slmi=newlmi
lmiterm([-Slmi 1 1 S],1,1)
lmiterm([Slmi 1 1 0],1)
lmisys=getlmis
其中:X 和 S 分别表示变量 X 和 ,而 BRL、Xpos 和 Slmi 则分别表示第 1、第 2 和第 3
个线性矩阵不等式。-Xpos 指的是第 2 个线性矩阵不等式的右边,-X 表示变量
S
X 的转
置。
lmiedit
线性矩阵不等式编辑器 lmiedit 是一个图形用户界面,它可以按符号方式直接确定线
性矩阵不等式系统。输入
lmiedit
出现一个具有一些可编辑文本区域和各种按钮的窗户。按以下步骤来确定一个线性矩阵不
Administrator
高亮
等式系统:
1.在文本区域的上半部分给出每一个矩阵变量的描述(名字和结构),其结构是通
过类型(S 表示对称块矩阵,R 表示无结构的长方矩阵,G 表示其他结构矩阵)和一个
“附加”的结构矩阵(类似于 lmivar 中的 struct)来刻画的。在文本编辑区,使用一行描
述一个变量。
2.在文本区的下半部分,按 MATLAB 的表示方式给出要描述的线性矩阵不等式。
例如,线性矩阵不等式
0<⎥⎦
⎤⎢⎣
⎡
−
+
IXB
XBXAXA
T
T
可以通过输入
[A’*X+X*A X*B;B’*X -1]<0
来描述。其中 X 是文本区上半部分描述矩阵变量 X 的变量名。一个线性矩阵不等式的描
述可能需要几行,但一行中最多只能描述一个线性矩阵不等式。
完成了线性矩阵不等式系统的描述后,可以通过按相应的按钮来完成以下的任务:
z 显示用于描述线性矩阵不等式的 lmivar/lmiterm 命令串(按钮 view commands);
反之,通过单击按钮 describe…可以将用一串 lmivar/lmiterm 命令定义的线性矩阵
不等式系统按 MATLAB 表示式显示。
z 将线性矩阵不等式的符号描述存为一个 MATLAB 语句串(按钮 save)。以后可以
通过按钮 load 重新显示这种描述。
z 可以从一个文件读一串 lmivar/lmiterm 命令(按钮 read),然后通过单击
“describe the matrix variables”或“describe the LMIs …”显示出由这些命令确定
的线性矩阵不等式系统的符号表示。
z 写一串用于描述一个特殊线性矩阵不等式系统的 lmivar/lmiterm 命令(按钮
write)。
z 通过按钮 create 产生线性矩阵不等式系统的内部表示,结果用一个以线性矩阵不
等式命名的 MATLAB 变量记录(如果线性矩阵不等式系统名是 mylmi,则其内部
表示用 MATLAB 变量 mylmi 记录)。内部表示 mylmi 可以被线性矩阵不等式求
解器或任何其他的线性矩阵不等式函数调用。
如同命令 lmiterm 一样,可以应用简捷的方法来输入线性矩阵不等式的表示式。例如
零块可以简单地输入 0,而不必定义其维数,类似地,单位矩阵只需输入数字 1 等。
lmiedit 尽管很一般,但它没有 lmiterm 灵活。以下是 lmiedit 的一些局限性:
z 在矩阵变量的两边不能使用括号。例如当 X 是一个变量名时,
(A*X+B)’*C+C’*(A*X+B)
是不允许的,而
(A+B)’*X+X’*(A+B)
则是可以的。
z 不允许出现循环和条件语句。
z 当把 lmiterm 命令转换成一个线性矩阵不等式的符号描述时,如果 lmiterm 的第 1
个分量不能确认就将出错。使用由 newlmi 和 lmivar 提供的线性矩阵不等式和变量
标识符可以避免这样的问题。
图 A.1 给出了用 lmiedit 描述例 1 中的线性矩阵不等式系统的窗口。
图 A.1 lmiedit 的图形界面
A.3 信息提取
线性矩阵不等式系统的完整描述是以一个叫做内部表示的向量储存在机器内的。LMI
工具箱提供了三个函数 lmiinfo、lminbr 和 matnbr,它们可以从内部表示向量中提取线性矩
阵不等式的相关信息,并以用户可读的方式显示出来。
lmiinfo
lmiinfo 是一种交互式工具,用以反映有关线性矩阵不等式系统的一些信息。这些信息
包括线性矩阵不等式的个数、矩阵变量的个数和它们的结构、每一个线性矩阵不等式块中
项的内容等。为了调用 lmiinfo,输入
lmiinfo(lmisys)
其中的 lmisys 是由 getlmis 产生的线性矩阵不等式系统的内部表示。
lminbr 和 matnbr
这两个函数给出了系统中线性矩阵不等式的个数和矩阵变量的个数。例如,为了得到
矩阵变量的个数,输入
matnbr(lmisys)
A.4 线性矩阵不等式求解器
LMI 工具箱提供了用于求解以下三个问题的线性矩阵不等式求解器(其中 表示决策
变量向量,即矩阵变量 中的独立变元构成的向量)。
x
kXX ,,1 "
z 可行性问题:
寻找一个 (或等价的:具有给定结构的矩阵 ),使得满足线性矩阵
不等式系统
NR∈x kXX ,,1 "
)()( xBxA <
相应的求解器是 feasp。
z 具有线性矩阵不等式约束的一个线性目标函数的最小化问题:
xc
x
Tmin
s.t. )()( xBxA <
相应的求解器是 mincx。
z 广义特征值的最小化问题:
λ
x
min
s.t. )()( xDxC < (9)
)(xB<0 (10)
)()( xBxA λ< (11)
相应的求解器是 gevp。
以下详细介绍这三个求解器的功能和使用方法。
feasp
求解器 feasp 的一般表达式如下:
[tmin,xfeas]=feasp(lmisys,options,target)
求解器 feasp 是通过求解如下的一个辅助凸优化问题
tmin
s.t. IxBxA t≤− )()(
来求解线性矩阵不等式系统 lmisys 的可行性问题。
这个凸优化问题的全局最优值用 tmin 表示,作为求解器 feasp 输出的第一个分量。如
果 tmin<0,则系统 lmisys 是可行的。当系统 lmisys 为可行时,求解器 feasp 输出的第二个
分量 xfeas 给出了该线性矩阵不等式系统决策变量的一个可行解。进而,应用 dec2mat 可
以得到系统 lmisys 矩阵变量的一个可行解。
求解器 feasp 的输入变量 target 为 tmin 设置了目标值,使得只要 tmin
R
2
1
2 Rx
N
i
i <∑
=
中,或者说向量 xfeas 的欧氏范数不超过 R。该参数的默认值是 。 910=R
可行域半径的设定可以避免产生具有很大数值的解 ,同时也可以加快计算过程,改
进数值稳定性。
x
z options(4):该参数用于加快迭代过程的结束,它提供了反映优化过程中迭代速度
和解的精度之间的一个折中指标。当该参数取值为一个正整数 J 时,表示在最后
的 次迭代中,如果每次迭代后 的减小幅度不超过 1%,则优化迭代过程就停
止。该参数的默认值是 10。
J t
z options(5):options(5)=1 表示不显示迭代过程中的数据,options(5)=0(默认值)则
相反。
将 options(i)设置为零相当于将相应的控制参数设置为默认值,也可以通过忽略该输
入变量来接受默认值。
例 3:求满足 IP > 的对称矩阵 P ,使得
0<+ 1T1 PAPA (12)
0<+ 2T2 PAPA (13)
0<+ 3T3 PAPA (14)
其中:
⎥⎦
⎤⎢⎣
⎡
−
−=⎥⎦
⎤⎢⎣
⎡
−
−=⎥⎦
⎤⎢⎣
⎡
−
−=
0.27.0
9.04.1
,
7.23.1
5.18.0
,
31
21
321 AAA
为了调用 feasp,我们首先确定线性矩阵不等式系统:
setlmis([])
P=lmivar(1,[2 1])
lmiterm([1 1 1 P],1,A1,’s’) % LMI #1
lmiterm([2 1 1 P],1,A2,’s’) % LMI #2
lmiterm([3 1 1 P],1,A3,’s’) % LMI #3
lmiterm([-4 1 1 P],1,1) % LMI #4: P
lmiterm([4 1 1 0],1) % LMI #4: I
lmis=getlmis
然后调用 feasp 来求该线性矩阵不等式系统的一个可行决策变量:
[tmin,xfeas]=feasp(lmis)
得到 tmin 。因此,线性矩阵不等式系统 lmis 是可行的。应用 dec2mat 1363.3−=
PP=dec2mat(lmis,xfeas,P)
得到问题的可行矩阵变量值:
⎥⎦
⎤⎢⎣
⎡=
1.1554.126
4.1268.270
P
在求解这个可行性问题的过程中,也可以附加一些约束,例如,要求矩阵 P 的
Frobenius 范数不超过 10,且 tmin 1−≤ 。可以通过调用
[tmin,xfeas]=feasp(lmis,[0,0,10,0,0],-1)
来达到这些附加要求。相应的结果是 tmin 1745.1−= ,相应的矩阵 P 的最大特征值是
。 6912.9)(max =Pλ
mincx
求解器 mincx 的一般表达式如下:
[copt,xopt]=mincx(lmisys,c,options,xinit,target)
问题中的线性矩阵不等式系统由 lmisys 表示,向量 c 和决策变量向量 x有相同的维
数。对于由矩阵变量表示的线性目标函数,可以应用函数 defcx 来得到适当的向量 c。函
数 mincx 返回到目标函数 的全局最优值 copt 和决策变量的最优解 xopt,相应的矩阵变
量的最优解可以应用函数 dec2mat 从 xopt 得到。
xcT
函数 mincx 的输入量中除了 lmisys 和 c 以外,其他的输入是可选择的。xinit 是最优解
xopt 的一个初始猜测(可以从矩阵变量 的给定值,通过使用 mat2dec 来导出
xinit)。当输入的 xinit 不是一个可行解时,它将被忽略;否则,则有可能加快问题求解的
过程。
kXX ,,1 "
target 是目标函数的一个设定目标,只要某个可行的 满足 target,求解过程就
停止。
x ≤xcT
options 是一个 5 维向量,它用来描述优化迭代过程中的一些控制参数:
z options(1):该参数确定了最优值 copt 所要求的精度(默认值是 )。 210−
z options(2):该参数设定优化迭代过程中允许的最大迭代次数(默认值是 100)。
z options(3):该参数设定了可行域的半径。有求解器 feasp 中的相应参数相同。
z options(4):该参数用于加快迭代过程的结束。当该参数取值为一个正整数 时,
表示在最后的 次迭代中。如果每次迭代后,目标函数 的减小幅度在给定的
精度内,则优化迭代过程就停止。该参数的默认值是 5。
J
J xcT
z options(5):options(5)=1 表示不显示迭代过程中的数据,options(5)=0(默认值)则
相反。
将 options(i)设置为零相当于将相应的控制参数设置为默认值,也可以通过忽略该输
入变量来接受默认值。
以下的例子说明了求解器 mincx 的使用方法。
例 4:考虑优化问题
)(Tracemin X
X
s.t. 0<+++ QXXBBXAXA TT
其中: X 是一个对称的矩阵变量,
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−−
−−−
−
=
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
=
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−−
−−
=
36120
1231
011
,
1
0
1
,
121
123
121
QBA
根据矩阵的 Schur 补性质,本例中的优化问题等价于
)(Tracemin X
X
s.t. 0<⎥⎦
⎤⎢⎣
⎡
−
++
IXB
XBQXAXA
T
T
由于 是)(Trace X X 的元的一个线性函数,因此以上的优化问题是一个具有线性矩阵
不等式约束的线性目标函数的最小化问题,从而可以应用求解器 mincx 来求解这个问题。
以下给出用 mincx 来求解该问题的过程。
1.定义线性矩阵不等式约束
setlmis([])
X=lmivar(1,[3 1]) % 变量 X,满对称的
lmiterm([1 1 1 X],1,A,’s’)
lmiterm([1 1 1 0],Q)
lmiterm([1 2 2 0],-1)
lmiterm([1 2 1 X],B’,1)
LMIs=getlmis
2.将目标函数 写成 ,其中)(Trace X xcT x是矩阵变量 X 中的独立元所构成的向量。
由于引进向量 的目的是要选择c X 的对角元,因此它可以作为相应于 IX = 的决策向量得
到,即
c=mat2dec(LMIs,eye(3))
事实上,函数 defcx 将提供一个确定这样的目标函数的更加系统化的方法。
3.调用 mincx 计算最小值 xopt,目标函数的全局最小值 copt=c’*xopt。
options=[1e-5,0,0,0,0]
[copt,xopt]=mincx(LMIs,c,options)
其中 1e-5 给定了所要求的关于 copt 的计算精度。
作为求解器 mincx 运行的结果,以下的信息将出现在屏幕上:
Solver for linear objective minimization under LMI constraints
Iterations : Best objective value so far
1
2 -8.511476
3 -13.063640
*** new lower bound: -34.023978
4 -15.768450
*** new lower bound: -25.005604
5 -17.123012
*** new lower bound: -21.306781
6 -17.882558
*** new lower bound: -19.819471
7 -18.339853
*** new lower bound: -19.189417
8 -18.552558
*** new lower bound: -18.919668
9 -18.646811
*** new lower bound: -18.803708
10 -18.687324
*** new lower bound: -18.753903
11 -18.705715
*** new lower bound: -18.732574
12 -18.712175
*** new lower bound: -18.723491
13 -18.714880
*** new lower bound: -18.719624
14 -18.716094
*** new lower bound: -18.717986
15 -18.716509
*** new lower bound: -18.717297
16 -18.716695
*** new lower bound: -18.716873
Result: feasible solution of required accuracy
best objective value: -18.716695
guaranteed relative accuracy: 9.50e-006
f-radius saturation: 0.000% of R = 1.00e+009
迭代的次数和当前这次迭代时 的最佳值分别在左列和右列中。注意,在第一次迭
代中没有一个对应的目标函数值,这表明满足约束条件的可行解 只是在第二次迭代时才
被找到。
xcT
x
4.mincx 也给出决策变量 xopt 的最优值。相应的矩阵变量最优值由以下函数给出
Xopt=dec2mat(LMIs,xopt,X)
由此可以得到:
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−
−−
−−
=
0771.62201.22046.2
2201.22855.65.8895
2046.28895.56.3542
optX
gevp
求解器 gevp 的一般表达式如下:
[lopt,xopt]=gevp(lmisys,nlfc,options,linit,xinit,target)
如果问题的线性矩阵不等式约束是可行的,则 gevp 给出了优化问题的全局最小值
lopt 和决策向量 x的最优解 xopt。相应的矩阵变量的最优解可以应用 dec2mat 得到。
输入分量 lmisys 表示当 1=λ 时由(6)~(8)构成的线性矩阵不等式系统。包含λ的
线性矩阵不等式系统称为线性分式约束。线性分式约束(8)的个数用 nlfc 表示。其他的
输入分量都是可选择的。
如果 是一个可行解,则通过令 linit=),( 00 xλ 0λ 、xinit= ,将 设置为 gevp 的
初始值。当 不是可行解时,这样的初始值设置不会被接受。target 的设定值表明了
只要当一个可行解
0x ),( 00 xλ
),( 00 xλ
),( xλ 满足 ≤λ target 时,迭代过程就停止。
可选择的输入量 options 是一个 5 维向量,它用来描述优化迭代过程中的一些控制参
数:
z options(1):该参数设定了最优值 lopt 所要求的精度(默认值是 )。 210−
z options(2):该参数设定优化迭代过程中允许的最大迭代次数(默认值是 100)。
z options(3):该参数设定了可行域的半径。与求解器 feasp 中的相应参数相同。
z options(4):该参数用于加快迭代过程的结束。当该参数取值为一个正整数 时,
表示在最后的
J
J 次迭代中,如果每次迭代后 λ的减小幅度在给定的精度内,则优
化迭代过程就停止。该参数的默认值是 5。
z options(5):options(5)=1 表示不显示迭代过程中的数据,options(5)=0(默认值)则
相反。
将 options(i)设置为零相当于将相应的控制参数设置为默认值,也可以通过忽略该输
入变量来接受默认值。
对广义特征值的最小化问题,在调用求解器 gevp 时,须遵循以下规则:
z 确定包含λ的线性矩阵不等式: )()( xBxA < (注意没有λ);
z 总是把 放在线性矩阵不等式系统的最后; )()( xBxA <
z 要求有约束 ,或者保证)(xB<0 )(xB<0 成立的任何其他约束。
根据以上要求,对于
0
00
0 >⎥⎦
⎤⎢⎣
⎡= )(,)()( 11 xBxBxB
的广义特征值优化问题,我们就不能直接应用求解器 gevp 来求解。为了克服这一困难,
可以通过引进矩阵变量Y ,并用
0
00
0 ><⎥⎦
⎤⎢⎣
⎡< )(),(,)( 1 xBxBYYxA λ
来代替
0>< )(),()( xBxBxA λ
而等价的新问题可以用 gevp 来求解。
例 5:对以下的 3 个系统
),3,2,1(),()( == itt i xAx�
其中矩阵 由例 3 给出。问题是寻找一个单一的 Lyapunov 函数
来验证给定的这 3 个系统的稳定性,同时使得衰减率
)3,2,1(, =iiA Pxxx T)( =V
t
V
d
)(d x− 最大化。这样一个问题等价
于如下的优化问题:
αmin
s.t. PI <
PPAPA α<+ 1T1
PPAPA α<+ 2T2
PPAPA α<+ 3T3
以下应用求解器 gevp 来求解该问题,为此,首先确定线性矩阵不等式系统:
setlmis([]);
P=lmivar(1,[2 1])
lmiterm([1 1 1 0],1) % P>I: I
lmiterm([-1 1 1 P],1,1) % P>I: P
lmiterm([2 1 1 P],1,A1,’s’) % LFC #1 (lhs)
lmiterm([-2 1 1 P],1,1) % LFC #1 (rhs)
lmiterm([3 1 1 P],1,A2,’s’) % LFC #2 (lhs)
lmiterm([-3 1 1 P],1,1) % LFC #2 (rhs)
lmiterm([4 1 1 P],1,A3,’s’) % LFC #3 (lhs)
lmiterm([-4 1 1 P],1,1) % LFC #3 (rhs)
lmis=getlmis
进而,通过以下命令来调用 gevp:
[alpha,popt]=gevp(lmis,3)
得到 alpha = -0.122 是该问题的最优值,因此相应的最大衰减率是 0.122,最优解是
⎥⎦
⎤⎢⎣
⎡
−
−=
64.1835.8
35.858.5
P
如何从决策变量到矩阵变量以及从矩阵变量到决策变量
当线性矩阵不等式由相应的矩阵变量描述时,线性矩阵不等式求解器涉及的是由这些
矩阵变量中的独立元所组成的决策向量 x。两个函数 mat2dec 和 dec2mat 可以实现这两种
变量之间的转换。
考虑一个具有三个矩阵变量 的线性矩阵不等式系统。给定这些变量的特
定值 X1、X2、X3,那么由 mat2dec 可以得到相应的决策向量的值:
321 XXX 、、
xdec=mat2dec(lmisys,X1,X2,X3)
如果 lmisys 后分量的个数和线性矩阵不等式系统 lmisys 中的矩阵变量个数不符,则系
统会提示一个出错信息。
这个函数在线性矩阵不等式求解器 mincx 或 gevp 的初始化中也是很有用的。例如,
给定 的一个初始猜测值,mat2dec 就形成了相应决策向量的初始值 xinit。 321 XXX 、、
反之,给定决策向量的一个值 xdec,那么可以通过函数 dec2mat 给出相应的第 k 个矩
阵的取值。例如,以下的表示式可以给出第 2 个矩阵变量的取值:
X2=dec2mat(lmisys,xdec,2)
函数 dec2mat 中的最后一个分量表明了要求的是第 2 个矩阵变量,这里也可以用
lmivar 定义的相应矩阵变量的变量名。
矩阵变量和决策变量的总数分别由 matnbr 和 decnbr 给出。另外,函数 decinfo 提供了
决策变量和矩阵变量之间关系的一些详细信息。
A.5 结果验证
LMI 工具箱提供了两个函数用于分析和验证一个线性矩阵不等式优化问题的结果。对
一个给定的决策向量的值,例如由线性矩阵不等式求解器给出的可行或最优解向量,函数
evallmi 求出线性矩阵不等式系统中所有变量项的值,进而应用 showlmi 给出特定线性矩阵
不等式的左边和右边。
在例 4 中,我们可以这样来验证由 mincx 得到的最优解 xopt 是否满足给定的约束条
件:
evlmi=evallmi(LMIs,xopt)
[lhs,rhs]=showlmi(evlmi,1)
第一个命令表示对给定的决策变量值 xopt,求取系统的值,第二个命令则显示第 1 个
线性矩阵不等式的左边和右边的矩阵值。这个不等式的成立与否可以通过
eig(lhs-rhs)
来检验,得到的结果是
ans=
-2.0387e-04
-3.9333e-05
-1.8917e-07
-4.6680e+01
由此可以看到 lhs-rhs 是负定的,因此 xopt 满足第 1 个线性矩阵不等式。
A.6 修改一个线性矩阵不等式系统
LMI 工具箱提供了函数 dellmi、delmvar 和 setmvar,它们可以用来修改一个已经确定
的线性矩阵不等式系统。以下来说明这些函数的用法。
Dellmi
用 dellmi 可以从一个线性矩阵不等式系统中删除一个完整的线性矩阵不等式。例如对
例 1 中的线性矩阵不等式系统,假定它由 lmisys 描述,现在要删除矩阵变量 X 的正定性
约束,这可以通过以下的命令来实现:
newsys=dellmi(lmisys,2)
其中第 2 个分量表明了要删除的是哪一个线性矩阵不等式。所导出的由两个线性矩阵不等
式构成的系统由 newsys 表示。
原来系统中的线性矩阵不等式标识符(由位置排列确定的)不因删除了其中的一些不
等式后而改变它们在新的线性矩阵不等式系统中的标识符。因此,在例 1 中删除了第 2 个
线性矩阵不等式后,
IS >
仍然是第 3 个线性矩阵不等式,尽管它在新的系统中已位于第 2 的位置。为了避免引起混
淆,更安全的方法是使用由 newlmi 给定的线性矩阵不等式名。如在例 1 中,BRL、Xpos
和 Slmi 依次是三个线性矩阵不等式的名称,则 Slmi 仍是通过
newsys=dellmi(lmisys,Xpos)
删除了第 2 个不等式后所得到的新系统中线性矩阵不等式 的名称。 IS >
delmvar
另一种修改一个线性矩阵不等式系统的方式是删除一个矩阵变量,即在所有包含该变
量的变量项中删除该矩阵变量。delmvar 提供了这样一种修改功能。例如,考虑线性矩阵
不等式系统
0<++++ IBWBWXAXA TTT
其中: ∈= TXX 44×R 和 是这个线性矩阵不等式中的两个矩阵变量。这个线性矩
阵不等式可以用以下的命令来描述:
42×∈RW
setlmis([])
X=lmivar(1,[4 1]) % X
W=lmivar(2,[2 4]) % S
lmiterm([1 1 1 X],1,A,’s’)
lmiterm([1 1 1 W],B,1,’s’)
lmiterm([1 1 1 0],1)
lmisys=getlmis
为了删除变量W ,输入命令
newsys=delmvar(lmisys,W)
得到的 newsys 描述了 Lyapunov 不等式
0IXAXA <++T
注意:delmvar 将自动删除所有的只依赖所删除矩阵变量的不等式。矩阵变量的变量名不
会因删除个别的矩阵变量而改变。
删除一个矩阵变量实际上等价于让这个矩阵变量等于零,这一功能也可以由下面的
setmvar 方便地实现。
setmvar
函数 setmvar 用于给某个矩阵变量赋值,一旦赋值,该变量就从原来的问题中消失
了,包含该变量的所有项都成为常数项。另外,这个函数在优化问题中也是很有用的,例
如可以通过固定某个变量,然后对其他变量进行优化。
考虑例 1,我们想要知道 的最大增益是否小