Solving feasibility problems results an error "Undefined function 'matchpairs' for input arguments of type 'double'."

[panqiincs]

I am designing a H infinity observer for linear discrete systems. Knowing the system matrices, I want to find an observer gain $L$ to satisfy the performance index. It is equivalent to the following feasibility problem:

I set a value for $\gamma$ and wrote following code to solve it:

%% Model description
n = 4;  % dimension of state x
l = 2;  % dimension of measured state y
m = 4;  % dimension of output z
p = 2;  % dimension of disturbance w
r = 1;  % dimension of input u

% n by n matrix
A = [0.9285, 0.3876,  0.0715,  0.0124;
    -0.3516, 0.9146,  0.3516,  0.0854;
     0.0071, 0.0012,  0.9929,  0.3988;
     0.0352, 0.0085, -0.0352, -0.9915];
% n by r matrix
B1 = [0.0013;
      0.0124;
      0.0799;
      0.3988];
% n by p matrix
B2 = [1, 0;
      0, 0;
      0, 1;
      0, 0];
% l by n matrix
C1 = [1, 0, 0, 0;
      0, 0, 1, 0];
% l by r matrix
D1 = [0; 0];
% l by p matrix
D2 = [1, 0; 0, 1];
% m by n matrix
C2 = eye(4);

%% Solve LMI problems
yalmip('clear')

gamma = 0.1;
P = sdpvar(n, n);
L = sdpvar(n, l, 'full');

M11 = (A+L*C1)'*P*(A+L*C1)-P;
M12 = (A+L*C1)'*P*(B2+L*D2);
M13 = C2';
M21 = M12';
M22 = (B2+L*D2)'*P*(B2+L*D2)-gamma*eye(p);
M23 = zeros(p, m);
M31 = M13';
M32 = M23';
M33 = -gamma*eye(m);
M = [M11, M12, M13;
     M21, M22, M23;
     M31, M32, M33];
Constraints = [P >= 0; M <= 0];

Objective = [];

options = sdpsettings('verbose', 1);

sol = optimize(Constraints, Objective, options);

if sol.problem == 0
    P = value(P);
    L = value(L);
    hinf_index = value(gamma);
else
    disp('Hmm, something went wrong!');
    sol.info
    yalmiperror(sol.problem)
end

It shows the following error:

>> hinf_observer_discrete
                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    0.000000e+00    4.000e+00    1.381e+00
Hmm, something went wrong!

ans =

    'Unknown problem (learn to debug) (Unknown problem (learn to debug) (Undefined function 'matchpairs' for input arguments of type 'double'.))'


ans =

    'Unknown problem (learn to debug) '

What's the problem and how to solve it? Many thanks!

[johanlofberg]

You are setting up a nonlinear nonconvex SDP (i.e. basically unsolvable today) and there are no solvers available for this practically. YALMIP has a built-in hack (called Slayer) where it tries to attack the problem using a nonlinear programming approach, and that algorithm uses the function matchpairs which was introduced in a rather recent version of MATLAB. Having said that, Slayer fails to solve the problem when running it here, not too surprisingly

[johanlofberg]

The built-in global solver BMIBNB (with the Slayer framework for the upper bound generation) manages to find a solution though

>> sol = optimize(Constraints, gamma, options);
* Starting YALMIP global branch & bound.
* Upper solver     : fmincon
* Lower solver     : MOSEK
* LP solver        : GUROBI
* -Extracting bounds from model
* -Perfoming root-node bound propagation
* -Calling upper solver (no solution found)
* -Branch-variables : 18
* -More root-node bound-propagation
* -Performing LP-based bound-propagation 
* -And some more root-node bound-propagation
* Warning: 14 branch variables are unbounded from below
* Warning: 17 branch variables are unbounded from above
* Starting the b&b process
 Node       Upper       Gap(%)       Lower     Open   Time
    1 :           Inf      NaN    0.00000E+00    2     1s    
    2 :           Inf      NaN    0.00000E+00    3     1s    
    3 :           Inf      NaN    0.00000E+00    4     2s    
    4 :           Inf      NaN    0.00000E+00    5     3s    
    5 :           Inf      NaN    0.00000E+00    6     4s    
    6 :   0.00000E+00     0.00    0.00000E+00    5     4s  Solution found by upper solver  
* Finished.  Cost: 0 (lower bound: 0, relative gap 0%)
* Termination with relative gap satisfied 
* Timing: 63% spent in upper solver (7 problems solved)
*         19% spent in lower solver (42 problems solved)
*         10% spent in LP-based domain reduction (180 problems solved)
*         1% spent in upper heuristics (11 candidates tried)

and by optimizing over gamma it can be pushed down significantly (but the lower bound does not appear to improve above the trivial bound 0 so problem is hard for the global solver, not too surprsingly as there are no strong variable bounds on the decision variables as warned about in the presolve

 Node       Upper       Gap(%)       Lower     Open   Time
 100 :   3.34924E-02     3.29    1.73469E-07   99   107s   

Although this model never can be made convex or linear (general output feedback) it can probably be improved (from the perspective of BMIBNB) as it currently is cubic, but I suspect some transformations at least can convert it to a quadratic model

[johanlofberg]

This comment

%% Solve LMI problems

Indicates that you don't understand what you are modelling. What does the L stand for...

[panqiincs]

Thanks for answering the question and pointing out my horriable mistake.