function [x, obj, Y, info] = runLMI(c, F, G, varargin) % RUNLMI Solve a semi-definite program given in dual (linear matrix % inequality) form % requires: SeDuMi % author: Toby Cubitt % license: GPL2 % % [x, OBJ, Y, INFO] = runLMI(c, F, G) solves the following SDP % (where >= denotes matrix inequality): % % minimize c.'*x % subject to sum_i F{k,i}*x(i) >= G{k} % % [x, OBJ, Y, INFO] = runLMI(c, F, G, A, b) solves the following SDP % with additional equality constraints: % % minimize c.'*x % subject to sum_i F{k,i}*x(i) >= G{k} % A*x == b % % The vector c defines the linear objective function, F is a % 2-dimensional cell array containing the LMI coefficient matrices, G is % a 1-dimensional cell array containing the LMI constant matrices, % matrix A and vector b specify the equality constraints. % % x is the optimal solution vector, OBJ is the value of the objective % function at that point, Y is the dual solution matrix. INFO is a % SeDuMi structure containing information about the outcome of the % numerical optimisation. % % [x, OBJ, Y, INFO] = runLMI(c, F, G, 'real') or % [x, OBJ, Y, INFO] = runLMI(c, F, G, A, b, 'real') assumes F and G are % real symmetric matrices. Otherwise, G and F can be Hermitian. % % If your problem only requires a single LMI, G is allowed to be a % single matrix rather than a cell array (F will be a % 1-dimensional cell array in this case). % % % See also RUNSDP, SEDUMI. % % Note that the optimisation problems solved by RUNLMI and RUNSDP are % not strictly dual to each other: if [x,f,Y] = runLMI(c,F,G), then % [Y,-f,-x] = runSDP(-G,F,c). %% Copyright (C) 2004-2011 Toby Cubitt %% %% This program is free software; you can redistribute it and/or %% modify it under the terms of the GNU General Public License %% as published by the Free Software Foundation; either version 2 %% of the License, or (at your option) any later version. %% %% This program is distributed in the hope that it will be useful, %% but WITHOUT ANY WARRANTY; without even the implied warranty of %% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the %% GNU General Public License for more details. %% %% You should have received a copy of the GNU General Public License %% along with this program; if not, write to the Free Software %% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, %% MA 02110-1301, USA. % sort out input arguments if isnumeric(G) if size(F,2) == 1 F = F'; end if size(F,1) ~= 1 error('inconsistent numbers of matrices in F and G'); end G = {G}; end R = 'complex'; if length(varargin) >= 3 R = varargin{3}; elseif length(varargin) >= 1 && ischar(varargin{1}) R = varargin{1}; end % get dimensions etc. nvars = length(c); nLMI = size(F,1); dim = zeros(nLMI,1); for k = 1:nLMI dim(k) = size(F{k,1},1); end %--- equality constraints --- At = sparse([]); bt = sparse([]); if length(varargin) >= 2 At = sparse(varargin{1}); bt = sparse(varargin{2}); K.f = size(At,1); end %--- LMI constraints --- K.s = dim; if ~ (strcmp(R,'real') || strcmp(R,'re')) K.scomplex = [1:nLMI]; end for k = 1:nLMI Atmp = sparse([]); for j = 1:nvars % !!! -1 because SeDuMi has LMI in form: G - x*F >= 0 !!! Atmp = [Atmp,vec(-sparse(F{k,j}))]; end At = [At;Atmp]; % !!! -1 because SeDuMi has LMI in form: G - x*F >= 0 !!! bt = [bt;vec(-sparse(G{k}))]; end %--- objective function --- % !!! -1 because SeDuMi maximizes dual rather than minimizing !!! c = -sparse(c); % ensure c is a column vector if size(c,1) == 1 c = c.'; end %% %--- display problem definition --- %% disp('objective function: c'); %% disp(full(c)); %% disp(''); %% A = []; %% if length(varargin) >= 2 %% A = varargin{1}; %% end %% disp('LMI:') %% disp('F'); %% for j = 1:nvars %% disp(full(mat(At(size(A,1)+1:size(A,1)+dim(1)^2,j)))); %% end %% for k = 2:nLMI %% for j = 1:nvars %% disp(''); %% disp(full(mat(At(size(A,1)+sum(dim(1:k-1).^2)+1:... %% size(A,1)+sum(dim(1:k).^2),j)))); %% end %% end %% disp(''); %% disp('G'); %% disp(full(mat(bt(size(A,1)+1:size(A,1)+dim(1)^2)))); %% for k = 2:nLMI %% disp(''); %% disp(full(mat(bt(sum(size(A,1)+dim(1:k-1).^2)+1:... %% size(A,1)+sum(dim(1:k).^2))))); %% end %% disp(''); %% if ~ (A == []) %% disp('equality constraint:') %% disp('A') %% disp(full(At(1:size(A,1),:))); %% disp('b'); %% disp(full(bt(1:size(A,1)))); %% end %--- run SDP --- [x,y,info] = sedumi(At, c, bt, K); Y{1} = mat(x(1:dim(1)^2)); for k = 2:nLMI Y{k} = mat(x(sum(dim(1:k-1).^2)+1:sum(dim(1:k).^2))); end x = y; obj = full(-c.'*x); % undo sign-change introduced into objective function