Decentralized linear controller synthesis via LMI: decLMI Class Definition

This class provides a tool to synthesize a stabilizing decentralized linear controller for discrete-time LTI systems. Functionalities are also given to achieve either robust or stochastic convergence to the origin.

Contents

Class description

This work is based on the paper "Synthesis of networked switching linear decentralized controllers" by D. Barcelli, D. Bernardini and A. Bemporad, 49th IEEE Conference on Decision and Control, 2010.

The inputs for the class are

All the methods solve an SDP problem and provide a feedback matrix gain K, and the matrix solution of the Riccati equation P. The SDP problems are dependent on the type of stability required. These methods are:

classdef decLMI

Output (read-only) properties

Mc is a structure which defines the Markov chain. It includes:

The other properties are the output of the SDP solution:

Each of the previous properties is a structure with fields:

    properties (SetAccess=private)
        Mc
        P
        K
        g
    end

Private properties

    properties (Access=private)
        Net
        A
        B
        Qx
        Qu
        X0
        xmax
        umax
        n
        m
        nLi
        Nconf
        conf
        tLambda_i
        tLambda
        li
        L
        ell
    end

    methods

Class constructor

obj=decLMI(Net,A,B,Qx,Qu,X0,xmax,umax,Mc)

Net : matrix of row size equal to number of actuators and columns size equal to number of states. Net(i,j) can be:

A : state matrix of the LTI system modeling the plant;

B : input matrix of the LTI system modeling the plant;

Qx : state weight matrix of the Riccati equation;

Qu : input weight matrix of the Riccati equation;

X0 : set of vertices of the polytope that defines the uncertainty of the initial state condition;

xmax : Euclidean norm constraint on state;

umax : Euclidean norm constraint on input;

Mc : two-states Markov chain that models the probability of losing a packet. Must be a structure with fields:

        function obj=decLMI(varargin)
            error(nargchk(8,9,nargin));

            Net=varargin{1};
            if 1%checkNet(Net)
                obj.Net=Net;
            else
                error('First argument must be a valid network topology');
            end
            clear Net;

            A=varargin{2};
            B=varargin{3};
            if (size(A,1)==size(A,2))&&(size(A,1)==size(B,1))
                obj.A=A;
                obj.B=B;
                [nx nu]=size(B);
            else
                error(['Second and Third arguments must be valid A and B '...
                    'state space matrix, respectively']);
            end
            clear A B;

            Qx=varargin{4};
            if (size(Qx,1)==size(Qx,2)) & (size(Qx,1)==nx)
                obj.Qx=Qx;
            else
                error(['Fourth argument must be a square matrix of the '...
                    'same size of A']);
            end
            clear Qx;

            Qu=varargin{5};
            if (size(Qu,1)==size(Qu,2)) & (size(Qu)==nu)
                obj.Qu=Qu;
            else
                error(['Fifth argument must be a square matrix of the '...
                    'same size of B']);
            end
            clear Qu;

            X0=varargin{6};
            if (size(X0,1)==nx)
                obj.X0=X0;
            else
                error(['Sixth argument must be an array of vertices of '...
                    'the polytope containing all initial conditions, '...
                    'hence each vertex must have same size as A']);
            end
            clear X0;

            xmax=varargin{7};
            if isnumeric(xmax)
                obj.xmax=xmax;
            else
                error(['Seventh argument is the norm of the states, '...
                    'hence must be a number']);
            end
            clear xmax;

            umax=varargin{8};
            if isnumeric(umax)
                obj.umax=umax;
            else
                error(['Eigth argument is the norm of the inputs, hence '...
                    'must be a number']);
            end
            clear umax;

            if (nargin==9)%check is in order for Mc
                obj.Mc=varargin{9};
            end

Computation of all possible Network configurations for all possible packet-loss configuration

            obj=obj.set_up_lambda();
            set(0,'RecursionLimit',1e4);
            tLambda_i = [];
            for i = 1:obj.m
                count=1;
                tLambda_i = define_tLambda_i(obj.Net,i,obj.nLi,tLambda_i,...
				obj.conf{i},count) ;
            end
            tLambda=[];
            tLambda = find_tilde_lambda(obj.Net,obj.conf,tLambda_i,tLambda);
            obj.tLambda=tLambda;
            obj.tLambda_i=tLambda_i;

Number of -1 in Net(i,:)

            obj.li = [];
            for i=1:obj.m
                obj.li(i) = sum(obj.Net(i,:)==-1);
            end

Total no. of -1 in Net

            obj.L = sum(obj.li);
            obj.ell = 2^obj.L;

            if (nargin==9)
                Mc=varargin{9};
                % transition matrix
                Mc.T = [Mc.q(1)   1-Mc.q(1);
                1-Mc.q(2) Mc.q(2)];
                obj.Mc=Mc;

                % create s0(i), s1(i)
                s0 = zeros(obj.ell,1);
                s1 = s0;
                for h=1:obj.ell
                    for i=1:obj.m
                        for j=1:obj.n
                            if obj.Net(i,j)==-1
                                s1(h) = s1(h) + obj.tLambda{h}(i,j);
                                s0(h) = s0(h) + (1-obj.tLambda{h}(i,j));
                            end
                        end
                    end
                end
                % emission matrix
                obj.Mc.E = zeros(2,obj.ell);
                for i=1:2
                    for j=1:obj.ell
                        obj.Mc.E(i,j)= obj.Mc.d(i)^(s0(j)) * ...
                            (1-obj.Mc.d(i))^(s1(j));
                    end
                end
            end
        end

Centralized LMI computation

Computes the solution of the SDP problem assuming the network to be fully connected and each link completely reliable. The goal is to give a reference for the performances that can be achieved via the decentralized methods.

        function obj = solve_centralized_lmi(obj)
            A=obj.A;
            B=obj.B;
            Qx=obj.Qx;
            Qu=obj.Qu;
            X0=obj.X0;
            xmax=obj.xmax;
            umax=obj.umax;

Formulate and solve an SDP for a centralized problem with no packet loss

            disp('Solving centralized SDP problem...');

            n = obj.n;
            m = obj.m;

Define SDP variables

            Q = sdpvar(n,n,'symmetric');
            Y = sdpvar(m,n);
            g = sdpvar(1,1);

Positive definiteness of Q

            F = set(Q>0);

Stability constraint

            F = F + set([Q          (A*Q+B*Y)'  (Qx^.5*Q)'  (Qu^.5*Y)'; ...
                         A*Q+B*Y    Q           zeros(n)    zeros(n,m); ...
                         Qx^.5*Q    zeros(n)    g*eye(n)    zeros(n,m); ...
                         Qu^.5*Y    zeros(m,n)  zeros(m,n)  g*eye(m)] >= 0);

Input constraints

            F = F + set([umax^2*eye(m)      Y;
                         Y'                 Q]>=0);

State constraints

            F = F + set([Q          (A*Q+B*Y)';
                         A*Q+B*Y    xmax^2*eye(n)]>=0);

            for i=1:size(X0,2)

Ellipsoid definition for every vertex of X0

                F = F + set([1          X0(:,i)'; ...
                            X0(:,i)     Q]>= 0);
            end

            sdpopt = sdpsettings('solver','sedumi');
            status_sol_cntr = solvesdp(F,g,sdpopt);

            if status_sol_cntr.problem~=0 && status_sol_cntr.problem~=4
                disp('\n\nCentralized problem infeasible\n\n')
                status_sol_cntr.problem
                return
            end

            disp('\nCentralized SDP problem solved!\n');

            invQ = inv(double(Q));
            g = double(g);
            P = g*invQ;
            K = double(Y)*invQ;

            OLeig = eig(A);
            CLeig = eig(A+B*K);

            figure
            r = rectangle('Position',[-1 -1 2 2],'Curvature',[1 1],...
                     'EdgeColor',[.5 .5 .5]);
            hold on
            plot(real(OLeig),imag(OLeig),'*b')
            plot(real(CLeig),imag(CLeig),'or')
            axis equal
            xlim = get(gca,'XLim');
            ylim = get(gca,'Ylim');
            plot([floor(xlim(1)) ceil(xlim(2))],[0 0],'Color',[.5 .5 .5]);
            plot([0 0],[-1.55 1.55],'Color',[.5 .5 .5]);
            title 'Centralized feedback: Eigenvalues'
            box

            obj.P.ci=P;
            obj.K.ci=K;
            obj.g.ci=g;
        end

Decentralized Ideal computation

Computes the solution considering only present links, but assuming that all of them are reliable. Basically does not account for miss-reception of packets containing mesurements.

        function obj = solve_dec_ideal_lmi(obj)
            A=obj.A;
            B=obj.B;
            Qx=obj.Qx;
            Qu=obj.Qu;
            X0=obj.X0;
            xmax=obj.xmax;
            umax=obj.umax;
            Net=obj.Net;

Formulate and solve an SDP for a decentralization with no packet loss

            disp('\nSolving decentralized SDP problem...\n');

            n = size(Net,2);
            m = size(Net,1);

Definition of zero components in K K(i,j) = 0 if state j is not available to compute input i

            K0 = Net;

Definition of the structure of Y

            Y0 = K0; % structure of Y

Definition of the variable Y

            Y = sdpvar(m,n);
            for i=1:m
                for j=1:n
                    if Y0(i,j)==0
                        Y(i,j) = 0;
                    end
                end
            end

Definition of the structure of Q

            Q0 = ones(n);
            for i=1:m
                for j=1:n
                    if K0(i,j)==0
                        for h=1:n
                            if K0(i,h)==1
                                Q0(h,j) = 0;
                                Q0(j,h) = 0;
                            end
                        end
                    end
                end
            end

Definition of the variable Q

            Q = sdpvar(n,n,'symmetric');
            for i=1:n
                for j=1:n
                    if Q0(i,j)==0
                        Q(i,j) = 0;
                    end
                end
            end
            g = sdpvar(1,1);

positive definiteness of Q

            F = set(Q>0);

Stability constraint

            F = F + set([Q  (A*Q+B*Y)'      (Qx^.5*Q)'  (Qu^.5*Y)'; ...
                A*Q+B*Y    Q               zeros(n)     zeros(n,m); ...
                Qx^.5*Q    zeros(n)        g*eye(n)     zeros(n,m); ...
                Qu^.5*Y    zeros(m,n)      zeros(m,n)   g*eye(m)] >= 0);

Input constraints

            F = F + set([umax^2*eye(m)      Y;
                Y'                 Q]>=0);

State constraints

            F = F + set([Q          (A*Q+B*Y)';
                A*Q+B*Y    xmax^2*eye(n)]>=0);

            for i=1:size(X0,2)

Ellipsoid definition for every vertex of X0

                F = F + set([1          X0(:,i)'; ...
                    X0(:,i)     Q]>= 0);
            end

            sdpopt = sdpsettings('solver','sedumi');
            status_sol_dec = solvesdp(F,g,sdpopt);

            if status_sol_dec.problem~=0 && status_sol_dec.problem~=4
                disp('\n\nDecentralized problem infeasible\n\n')
                return
            end

            disp('\nDecentralized SDP problem solved!\n');

            invQ = inv(double(Q));
            g = double(g);
            P = g*invQ;
            K = double(Y)*invQ;

            OLeig = eig(A);
            CLeig = eig(A+B*K);

            figure
            r = rectangle('Position',[-1 -1 2 2],'Curvature',[1 1],...
                     'EdgeColor',[.5 .5 .5]);
            hold on
            plot(real(OLeig),imag(OLeig),'*b')
            plot(real(CLeig),imag(CLeig),'or')
            axis equal
            xlim = get(gca,'XLim');
            ylim = get(gca,'Ylim');
            plot([floor(xlim(1)) ceil(xlim(2))],[0 0],'Color',[.5 .5 .5]);
            plot([0 0],[-1.55 1.55],'Color',[.5 .5 .5]);
            title 'Decentralized ideal feedback: Eigenvalues'
            box

            obj.P.di=P;
            obj.K.di=K;
            obj.g.di=g;
        end

Decentralized Robust computation

Computes a decentralized controller which provides stability and constraints satisfaction for any possible occurrence of the packet dropouts.

        function obj = solve_dec_lossy_lmi(obj)
            A=obj.A;
            B=obj.B;
            Qx=obj.Qx;
            Qu=obj.Qu;
            X0=obj.X0;
            xmax=obj.xmax;
            umax=obj.umax;
            Net=obj.Net;

Formulate and solve an SDP for a decentralization with packet-loss

            disp('\nSolving dec. lossy SDP problem...\n');
            n = size(Net,2);
            m = size(Net,1);

In case there is no wireless links in a row, a fixed structure have to be applied: K0fix. It should not be necessary to allocate all K0fix, but int his way it is useful to manage indices

            K0fix = sdpvar(m,n);
            nLi = obj.nLi;
            for i=1:m
                Nconf(i) = 2^(length(nLi{i}));
                if Nconf(i)>1
                    conf{i} = dec2bin(Nconf(i)-1);
                else

There is no -1 in this line

                    conf{i} = '-';

Set cooherently line in Kofix

                    for j = 1:n
                        if Net(i,j)==0
                            K0fix(i,j) = 0;
                        end
                    end
                end

            end

Define a memory structure that contains all possible configurations for each line of Net and its implications in K0{i,j}, where i means the i-th line of K0 and j means one of the 2^nLi{i} possible configurations

            K0i = [];
            set(0,'RecursionLimit',1e4);

In order to exploit the exponetial structure, the set of all possible configurations is generated using some recursive function that are not methods of the class.

            for i = 1:m
                count=1;
                K0i = define_K0_i(Net,i,nLi,K0i,conf{i},count) ;
            end
            Y = [];

            Y = def_lmis_K0i(K0i,nLi,K0fix,Y,conf,Net);

            Q = sdpvar(n,n,'symmetric');
            for w=1:length(Y)
                for i=1:m
                    for j=1:n
                        if double(Y{w}(i,j))==0
                            for h=1:n
                                if isnan(double(Y{w}(i,h)))
                                    Q(h,j)=0;
                                    Q(j,h)=0;
                                end
                            end
                        end
                    end
                end
            end
            g = sdpvar(1,1);
            NN = length(Y);

Positive definiteness of Q

            F = set(Q > 0);
            for j=1:size(X0,2)

Ellipsoid definition for every vertex of X0

                F = F + set([1          X0(:,j)'; ...
                    X0(:,j)     Q       ]>= 0);
            end

For all posssible configurations

            for i=1:NN

Stability constraint

                F = F + set([Q  (A*Q+B*Y{i})' (Qx^.5*Q)'   (Qu^.5*Y{i})'; ...
                    A*Q+B*Y{i} Q             zeros(n)     zeros(n,m); ...
                    Qx^.5*Q    zeros(n)      g*eye(n)     zeros(n,m); ...
                    Qu^.5*Y{i} zeros(m,n)    zeros(m,n)   g*eye(m)] >= 0);

Input constraints

                F = F + set([umax^2*eye(m)     Y{i};
                    Y{i}'              Q  ]>=0);

State constraints

                F = F + set([Q             (A*Q+B*Y{i})';
                    A*Q+B*Y{i}    xmax^2*eye(n)]>=0);
            end
            sdpopt = sdpsettings('solver','sedumi');
            sol = solvesdp(F,g,sdpopt);
            if sol.problem~=0 && sol.problem~=4
                disp('\n\nDecentralized lossy problem infeasible\n\n')
                sol.problem
                return
            end
            g = double(g);
            invQ = inv(double(Q));
            P = g*invQ;
            for i=1:NN
                K{i} = double(Y{i})*invQ;
            end
            for i=1:NN
                disp(['Closed loop max eig: ' ...
                    num2str(max(abs(eig(A+B*K{i}))))]);
            end
            OLeig = eig(A);
            CLeig=[];
            for i=1:NN
                CLeig = [CLeig;eig(A+B*K{i})];
            end

            figure
            r = rectangle('Position',[-1 -1 2 2],'Curvature',[1 1],...
                'EdgeColor',[.5 .5 .5]);
            hold on
            plot(real(OLeig),imag(OLeig),'*b')
            plot(real(CLeig),imag(CLeig),'or')
            axis equal
            xlim = get(gca,'XLim');
            ylim = get(gca,'Ylim');
            plot([floor(xlim(1)) ceil(xlim(2))],[0 0],'Color',[.5 .5 .5]);
            plot([0 0],[-1.55 1.55],'Color',[.5 .5 .5]);
            title 'Decentralized lossy feedback: Eigenvalues'
            box

            obj.P.dl=P;
            obj.K.dl=K;
            obj.g.dl=g;
        end

Decentralized Stochastic computation

Computes a decentralized controller that exploits the available knowledge on the dropouts probabiliy distributions. A two-states Markov chain is used to model the packet losses. Mean-square stability is guaranteed.

        function obj = solve_dec_stoch_lmi(obj)%v4
            if ~isempty(obj.Mc)
                Net=obj.Net;
                A=obj.A;
                B=obj.B;
                Qx=obj.Qx;
                Qu=obj.Qu;
                X0=obj.X0;
                xmax=obj.xmax;
                umax=obj.umax;
                Mc=obj.Mc;

Formulate and solve an SDP for a decentralization with packet loss and convergence in mean-square

                disp('\nSolving decentralized stoch. SDP problem...\n');
                n = size(Net,2);
                m = size(Net,1);

In case there are no wireless links in a row, a fixed structure have to be applied: K0fix.

                K0fix1 = sdpvar(m,n);
                K0fix2 = sdpvar(m,n);
                obj = obj.find_all_lossy();
                nLi=obj.nLi;
                for i=1:m
                    Nconf(i) = 2^(length(nLi{i}));
                    if Nconf(i)>1
                        conf{i} = dec2bin(Nconf(i)-1);
                    else

There is no -1 in this line

                        conf{i} = '-';

Set cooherently Kofix's line

                        for j = 1:n
                            if Net(i,j)==0
                                K0fix1(i,j) = 0;
                                K0fix2(i,j) = 0;
                            end
                        end
                    end

                end

Define a memory structure that contains all possible configurations for each line of Net and its implications in K0{i,j}, where i means the i-th line of K0 and j means one of the 2^nLi{i} possible configurations

                K0i = [];
                set(0,'RecursionLimit',1e4);
                for i = 1:m
                    count=1;
                    K0i = define_K0_i(Net,i,nLi,K0i,conf{i},count);
                end
                Y = [];
                Y1= [];
                Y2= [];
                Q=[];
                Y1 = def_lmis_K0i(K0i,nLi,K0fix1,Y1,conf,Net);

For all possible states of the Markov chain

                for ww=1:2
                    Q{ww} = sdpvar(n,n,'symmetric');

For all possible network configurations in that Markov chain state

                    for w=1:length(Y1)
                        for i=1:m
                            for j=1:n
                                if double(Y1{w}(i,j))==0
                                    for h=1:n
                                        if isnan(double(Y1{w}(i,h)))
                                            Q{ww}(h,j)=0;
                                            Q{ww}(j,h)=0;
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
                Y2 = def_lmis_K0i(K0i,nLi,K0fix2,Y2,conf,Net);
                Y = {Y1' Y2'};

                l = length(Y1);

Total number of lossy links

                L=log2(l);

                conf=[];
                for i=1:L
                    conf=[conf '1'];
                end

                d=Mc.d;

                e=Mc.E;

                re=e;
                for i=1:2
                    for j=1:l
                        re(i,j)=sqrt(e(i,j));
                    end
                end


                g = sdpvar(1,1);

                rhox = sdpvar(1,1);
                rhou = sdpvar(1,1);

Positive definiteness of Q in state 1

                F = set(Q{1} > 0);

Positive definiteness of Q in state 2

                F = set(Q{2} > 0);

                F = F + set(g > 0);
                F = F + set(rhox > 0);
                F = F + set(rhou > 0);

                Qrho=1e3*[1 1];


                for jj=1:size(X0,2)

Ellipsoid definition for every vertex of X0

                    F = F + set([1          X0(:,jj)'; ...
                        X0(:,jj)    Q{1}     ]>=0);

Ellipsoid definition for every vertex of X0

                    F = F + set([1          X0(:,jj)'; ...
                        X0(:,jj)    Q{2}     ]>=0);
                end

                for j=1:2

Define the constraint as [Qj [CC' CC1'];[CC;CC1] DD]

                    CC=[];
                    for i=1:l
                        CC=[CC ; sqrt(Mc.T(j,1))*re(j,i)*(A*Q{j}+B*Y{j}{i})];
                    end
                    for i=1:l
                        CC=[CC ; sqrt(Mc.T(j,2))*re(j,i)*(A*Q{j}+B*Y{j}{i})];
                    end
                    CC=[CC ; Qx^.5*Q{j}];
                    for i=1:l
                        CC=[CC ; re(j,i)*Qu^.5*Y{j}{i}];
                    end

                    str=[];
                    for i=1:l
                        str=[str 'Q{1},'];
                    end
                    for i=1:l
                        str=[str 'Q{2},'];
                    end
                    str=[str 'g*eye(n),'];
                    for i=1:l
                        str=[str 'g*eye(m),'];
                    end
                    str(end)=')';
                    eval(['DD=blkdiag(' str ';']);

                    F = F + set([Q{j} CC';CC DD]>=0);
                end

Implement hard constraints

                F = F + set(rhox==0);
                F = F + set(rhou==0);

                sdpopt = sdpsettings('solver','sedumi');
                status_sol_dec = solvesdp(F,g+Qrho*[rhox;rhou],sdpopt);

                if status_sol_dec.problem~=0 && status_sol_dec.problem~=4
                    disp('\n\nDecentralized stoch. problem infeasible\n\n')
                    return
                end

                disp('\nDecentralized stoch. SDP problem solved!\n');

                g = double(g);
                for j=1:2
                    invQ = inv(double(Q{j}));
                    P{j} = g*invQ;
                    for i=1:l
                        K{i,j} = double(Y{j}{i})*invQ;
                    end
                end


                OLeig = eig(A);
                CLeig=[];
                for j=1:2
                    for i=1:l
                        CLeig = [CLeig;eig(A+B*K{i,j})];
                    end
                end

                figure
                r = rectangle('Position',[-1 -1 2 2],'Curvature',[1 1],...
                    'EdgeColor',[.5 .5 .5]);
                hold on
                plot(real(OLeig),imag(OLeig),'*b')
                plot(real(CLeig),imag(CLeig),'or')
                axis equal
                xlim = get(gca,'XLim');
                ylim = get(gca,'Ylim');
                plot([floor(xlim(1)) ceil(xlim(2))],[0 0],'Color',[.5 .5 .5]);
                plot([0 0],[-1.55 1.55],'Color',[.5 .5 .5]);
                title 'Decentralized stochastic feedback: Eigenvalues'
                box

                obj.P.ds=P;
                obj.K.ds=K;
                obj.g.ds=g;
            else
                disp('Markov chain is needed. Computation is not performed.')
            end
        end
    end

Private methods

Not described below since their use is depicted in the code settion where the method is called.

    methods (Access=private)

        function obj=set_up_lambda(obj)
            obj.n = size(obj.Net,2);
            obj.m = size(obj.Net,1);
            obj = obj.find_all_lossy();
            for i=1:obj.m
                obj.Nconf(i) = 2^(length(obj.nLi{i}));
                if obj.Nconf(i)>1
                    obj.conf{i} = dec2bin(obj.Nconf(i)-1);
                else
                    obj.conf{i} = '-'; % there is no -1 in this line
                end
            end
        end

        function obj = find_all_lossy(obj)
            obj.nLi=[];
            for i=1:size(obj.Net,1)
                obj.nLi{i} = find(obj.Net(i,:)==-1);
            end
        end

     end
end