Documentation |
This last section gives a few hints for making the most out of the LMI Lab. It is directed toward users who are comfortable with the basics, as described in Tools for Specifying and Solving LMIs.
Fairly complex matrix variable structures and interdependencies can be specified with lmivar. Recall that the symmetric block-diagonal or rectangular structures are covered by Types 1 and 2 of lmivar provided that the matrix variables are independent. To describe more complex structures or correlations between variables, you must use Type 3 and specify each entry of the matrix variables directly in terms of the free scalar variables of the problem (the so-called decision variables).
With Type 3, each entry is specified as either 0 or ±x_{n} where x_{n} is the n-th decision variable. The following examples illustrate how to specify nontrivial matrix variable structures with lmivar. First, consider the case of uncorrelated matrix variables.
Suppose that the problem variables include a 3-by-3 symmetric matrix X and a 3-by-3 symmetric Toeplitz matrix:
$$Y=\left(\begin{array}{ccc}{y}_{1}& {y}_{2}& {y}_{3}\\ {y}_{2}& {y}_{1}& {y}_{2}\\ {y}_{3}& {y}_{2}& {y}_{1}\end{array}\right).$$
The variable Y has three independent entries, hence involves three decision variables. Since Y is independent of X, these decision variables should be labeled n + 1, n + 2, n + 3 where n is the number of decision variables involved in X. To retrieve this number, define the variable X (Type 1) by
setlmis([]) [X,n] = lmivar(1,[3 1])
The second output argument n gives the total number of decision variables used so far (here n = 6). Given this number, Y can be defined by
Y = lmivar(3,n+[1 2 3;2 1 2;3 2 1])
or equivalently by
Y = lmivar(3,toeplitz(n+[1 2 3]))
where toeplitz is a standard MATLAB^{®} function. For verification purposes, we can visualize the decision variable distributions in X and Y with decinfo:
lmis = getlmis decinfo(lmis,X) ans = 1 2 4 2 3 5 4 5 6 decinfo(lmis,Y) ans = 7 8 9 8 7 8 9 8 7
The next example is a problem with interdependent matrix variables.
Consider three matrix variables X, Y, Z with structure
$$X=\left(\begin{array}{cc}x& 0\\ 0& y\end{array}\right),\text{\hspace{1em}}Y=\left(\begin{array}{cc}z& 0\\ 0& t\end{array}\right),\text{\hspace{1em}}Z=\left(\begin{array}{cc}0& -x\\ -t& 0\end{array}\right)$$
where x, y, z, t are independent scalar variables. To specify such a triple, first define the two independent variables X and Y (both of Type 1) as follows:
setlmis([]) [X,n,sX] = lmivar(1,[1 0;1 0]) [Y,n,sY] = lmivar(1,[1 0;1 0])
The third output of lmivar gives the entry-wise dependence of X and Y on the decision variables (x_{1}, x_{2}, x_{3}, x_{4}) := (x, y, z, t):
sX = 1 0 0 2 sY = 3 0 0 4
Using Type 3 of lmivar, you can now specify the structure of Z in terms of the decision variables x_{1} = x and x_{4} = t:
[Z,n,sZ] = lmivar(3,[0 -sX(1,1);-sY(2,2) 0])
Since sX(1,1) refers to x_{1} while sY(2,2) refers to x_{4}, this defines the variable
$$Z=\left(\begin{array}{cc}0& -{x}_{1}\\ -{x}_{4}& 0\end{array}\right)=\left(\begin{array}{cc}0& -x\\ -t& 0\end{array}\right)$$
as confirmed by checking its entry-wise dependence on the decision variables:
sZ = 0 -1 -4 0
The LMI solvers are written for real-valued matrices and cannot directly handle LMI problems involving complex-valued matrices. However, complex-valued LMIs can be turned into real-valued LMIs by observing that a complex Hermitian matrix L(x) satisfies
L(x) < 0
if and only if
$$\left(\begin{array}{cc}\mathrm{Re}\left(L\left(x\right)\right)& \mathrm{Im}\left(L\left(x\right)\right)\\ -\mathrm{Im}\left(L\left(x\right)\right)& \mathrm{Re}\left(L\left(x\right)\right)\end{array}\right)<0.$$
This suggests the following systematic procedure for turning complex LMIs into real ones:
Decompose every complex matrix variable X as
X = X_{1} + jX_{2}
where X_{1} and X_{2} are real
Decompose every complex matrix coefficient A as
A = A_{1} + jA_{2}
where A_{1} and A_{2} are real
Carry out all complex matrix products. This yields affine expressions in X_{1}, X_{2} for the real and imaginary parts of each LMI, and an equivalent real-valued LMI is readily derived from the above observation.
For LMIs without outer factor, a streamlined version of this procedure consists of replacing any occurrence of the matrix variable X = X_{1} + jX_{2 }by
$$\left(\begin{array}{cc}{X}_{1}& {X}_{2}\\ -{X}_{2}& {X}_{1}\end{array}\right)$$
and any fixed matrix A = A_{1} + jA_{2}, including real ones, by
$$\left(\begin{array}{cc}{A}_{1}& {A}_{2}\\ -{A}_{2}& {A}_{1}\end{array}\right).$$
For instance, the real counterpart of the LMI system
M^{H}XM < X, X = X^{H} > I | (4-11) |
reads (given the decompositions M = M_{1} + jM_{2} and X = X_{1} + jX_{2} with Mj, Xj real):
$$\begin{array}{c}{\left(\begin{array}{cc}{M}_{1}& {M}_{2}\\ -{M}_{2}& {M}_{1}\end{array}\right)}^{T}\left(\begin{array}{cc}{X}_{1}& {X}_{2}\\ -{X}_{2}& {X}_{1}\end{array}\right)\left(\begin{array}{cc}{M}_{1}& {M}_{2}\\ -{M}_{2}& {M}_{1}\end{array}\right)<\left(\begin{array}{cc}{X}_{1}& {X}_{2}\\ -{X}_{2}& {X}_{1}\end{array}\right)\\ \left(\begin{array}{cc}{X}_{1}& {X}_{2}\\ -{X}_{2}& {X}_{1}\end{array}\right)<I.\end{array}$$
Note that X = X^{H} in turn requires that $${X}_{1}={X}_{1}^{H}$$ and $${X}_{2}+{X}_{2}^{T}=0$$. Consequently, X_{1} and X_{2} should be declared as symmetric and skew- symmetric matrix variables, respectively.
Assuming, for instance, that M ∊ C^{5}^{×}^{5}, the LMI system (Equation 4-11) would be specified as follows:
M1=real(M), M2=imag(M) bigM=[M1 M2;-M2 M1] setlmis([]) % declare bigX=[X1 X2;-X2 X1] with X1=X1' and X2+X2'=0: [X1,n1,sX1] = lmivar(1,[5 1]) [X2,n2,sX2] = lmivar(3,skewdec(5,n1)) bigX = lmivar(3,[sX1 sX2;-sX2 sX1]) % describe the real counterpart of (1.12): lmiterm([1 1 1 0],1) lmiterm([-1 1 1 bigX],1,1) lmiterm([2 1 1 bigX],bigM',bigM) lmiterm([-2 1 1 bigX],1,1) lmis = getlmis
Note the three-step declaration of the structured matrix variable bigX,
$$bigX=\left(\begin{array}{cc}{X}_{1}& {X}_{2}\\ -{X}_{2}& {X}_{1}\end{array}\right).$$
Specify X_{1} as a (real) symmetric matrix variable and save its structure description sX1 as well as the number n1 of decision variables used in X_{1}.
Specify X_{2} as a skew-symmetric matrix variable using Type 3 of lmivar and the utility skewdec. The command skewdec(5,n1) creates a 5-by–5 skew-symmetric structure depending on the decision variables n1 + 1, n1 + 2,...
Define the structure of bigX in terms of the structures sX1 and sX2 of X_{1} and X_{2}.
See Structured Matrix Variables for more details on such structure manipulations.
The LMI solver mincx minimizes linear objectives of the form c^{T}x where x is the vector of decision variables. In most control problems, however, such objectives are expressed in terms of the matrix variables rather than of x. Examples include Trace(X) where X is a symmetric matrix variable, or u^{T}Xu where u is a given vector.
The function defcx facilitates the derivation of the c vector when the objective is an affine function of the matrix variables. For the sake of illustration, consider the linear objective
$$\text{Trace}\left(X\right)+{x}_{0}^{T}P{x}_{0}$$
where X and P are two symmetric variables and x_{0} is a given vector. If lmsisys is the internal representation of the LMI system and if x_{0}, X, P have been declared by
x0 = [1;1] setlmis([]) X = lmivar(1,[3 0]) P = lmivar(1,[2 1]) : : lmisys = getlmis
the c vector such that $${c}^{T}x=\text{Trace}\left(X\right)+{x}_{0}^{T}P{x}_{0}$$ can be computed as follows:
n = decnbr(lmisys) c = zeros(n,1) for j=1:n, [Xj,Pj] = defcx(lmisys,j,X,P) c(j) = trace(Xj) + x0'*Pj*x0 end
The first command returns the number of decision variables in the problem and the second command dimensions c accordingly. Then the for loop performs the following operations:
Evaluate the matrix variables X and P when all entries of the decision vector x are set to zero except xj: = 1. This operation is performed by the function defcx. Apart from lmisys and j, the inputs of defcx are the identifiers X and P of the variables involved in the objective, and the outputs Xj and Pj are the corresponding values.
Evaluate the objective expression for X:= Xj and P:= Pj. This yields the j-th entry of c by definition.
In our example the result is
c = 3 1 2 1
Other objectives are handled similarly by editing the following generic skeleton:
n = decnbr( LMI system ) c = zeros(n,1) for j=1:n, [ matrix values ] = defcx( LMI system,j, matrix identifiers) c(j) = objective(matrix values) end
When solving LMI problems with feasp, mincx, or gevp, it is possible to constrain the solution x to lie in the ball
x^{T}x < R^{2}
where R > 0 is called the feasibility radius. This specifies a maximum (Euclidean norm) magnitude for x and avoids getting solutions of very large norm. This may also speed up computations and improve numerical stability. Finally, the feasibility radius bound regularizes problems with redundant variable sets. In rough terms, the set of scalar variables is redundant when an equivalent problem could be formulated with a smaller number of variables.
The feasibility radius R is set by the third entry of the options vector of the LMI solvers. Its default value is R = 109. Setting R to a negative value means "no rigid bound," in which case the feasibility radius is increased during the optimization if necessary. This "flexible bound" mode may yield solutions of large norms.
The LMI solvers used in the LMI Lab are based on interior-point optimization techniques. To compute feasible solutions, such techniques require that the system of LMI constraints be strictly feasible, that is, the feasible set has a nonempty interior. As a result, these solvers may encounter difficulty when the LMI constraints are feasible but not strictly feasible, that is, when the LMI
L(x) ≤ 0
has solutions while
L(x) < 0
has no solution.
For feasibility problems, this difficulty is automatically circumvented by feasp, which reformulates the problem:
Find x such that
L(x) ≤ 0 | (4-12) |
as:
Minimize t subject to
Lx < t × I. | (4-13) |
In this modified problem, the LMI constraint is always strictly feasible in x, t and the original LMI Equation 4-12 is feasible if and only if the global minimum t_{min} of Equation 4-12 satisfies
t_{min} ≤ 0
For feasible but not strictly feasible problems, however, the computational effort is typically higher as feasp strives to approach the global optimum t_{min} = 0 to a high accuracy.
For the LMI problems addressed by mincx and gevp, nonstrict feasibility generally causes the solvers to fail and to return an "infeasibility" diagnosis. Although there is no universal remedy for this difficulty, it is sometimes possible to eliminate underlying algebraic constraints to obtain a strictly feasible problem with fewer variables.
Another issue has to do with homogeneous feasibility problems such as
A^{T}P + P A < 0, P > 0
While this problem is technically well-posed, the LMI optimization is likely to produce solutions close to zero (the trivial solution of the nonstrict problem). To compute a nontrivial Lyapunov matrix and easily differentiate between feasibility and infeasibility, replace the constraint P > 0-by-P > αI with α > 0. Note that this does not alter the problem due to its homogeneous nature.
Consider the generalized eigenvalue minimization problem
Minimize λ subject to
A(x) < λB(x), B(x) > 0, C(x) <0. | (4-14) |
Technically, the positivity of B(x) for some x ∊ R^{n} is required for the well-posedness of the problem and the applicability of polynomial-time interior-point methods. Hence problems where
$$B\left(x\right)=\left(\begin{array}{cc}{B}_{1}\left(x\right)& 0\\ 0& 0\end{array}\right)$$
with B_{1}(x) > 0 strictly feasible, cannot be directly solved with gevp. A simple remedy consists of replacing the constraints
A(x) < B(x), B(x) > 0
by
$$A\left(x\right)<\left(\begin{array}{cc}Y& 0\\ 0& 0\end{array}\right),\text{\hspace{1em}}Y<\lambda {B}_{1}\left(x\right),\text{\hspace{1em}}{B}_{1}\left(x\right)>0$$
where Y is an additional symmetric variable of proper dimensions. The resulting problem is equivalent to Equation 4-14 and can be solved directly with gevp.
As explained in Tools for Specifying and Solving LMIs, the term-oriented description of LMIs used in the LMI Lab typically leads to higher efficiency than the canonical representation
A_{0} + x_{1}A_{1} + ... + x_{N}A_{N} < 0. | (4-15) |
This is no longer true, however, when the number of variable terms is nearly equal to or greater than the number N of decision variables in the problem. If your LMI problem has few free scalar variables but many terms in each LMI, it is therefore preferable to rewrite it as Equation 4-15 and to specify it in this form. Each scalar variable x_{j} is then declared independently and the LMI terms are of the form x_{j}A_{j}.
If M denotes the total row size of the LMI system and N the total number of scalar decision variables, the flop count per iteration for the feasp and mincx solvers is proportional to
N^{3} when the least-squares problem is solved via. Cholesly factorization of the Hessian matrix (default) [2]
M-by-N^{2} when numerical instabilities warrant the use of QR factorization instead
While the theory guarantees a worst-case iteration count proportional to M, the number of iterations actually performed grows slowly with M in most problems. Finally, while feasp and mincx are comparable in complexity, gevp typically demands more computational effort. Make sure that your LMI problem cannot be solved with mincx before using gevp.
In many output-feedback synthesis problems, the design can be performed in two steps:
Compute a closed-loop Lyapunov function via LMI optimization.
Given this Lyapunov function, derive the controller state-space matrices by solving an LMI of the form
M + P^{T}XQ + Q^{T}X^{T}P < 0 | (4-16) |
where M, P, Q are given matrices and X is an unstructured m-by-n matrix variable.
It turns out that a particular solution Xc of Equation 4-16 can be computed via simple linear algebra manipulations [1]. Typically, Xc corresponds to the center of the ellipsoid of matrices defined by Equation 4-16.
The function basiclmi returns the "explicit" solution Xc:
Xc = basiclmi(M,P,Q)
Since this central solution sometimes has large norm, basiclmi also offers the option of computing an approximate least-norm solution of Equation 4-16. This is done by
X = basiclmi(M,P,Q,'Xmin')
and involves LMI optimization to minimize ∥X ∥.