Need help in Proper modeling to find maximum ellipsoid inscribed in a polytope

[Sunandan Adhikary]

Here is my sample code to find the maximum inscribed ellipsoid: {x : (x-c)'*P^(-1)*(x-c)<=1}

contained in a polytope: {x : A*x <= b}, using YALMIP.

A = [ eye(3); -eye(3) ];

b = ones(6,1);

yref = zeros(3,1);

y = sdpvar(3, 1);

[P, ~] = maxInscribedEllipsoid(A, b, yref);

ellipsoid_constr = (y-c)*P*(y-c) == 1;

figure();

constr_poly = A*y - b <= 0;

confg = sdpsettings('verbose', 0);

confg.plot.linewidth = 2;

confg.plot.shade = 0.2;

plot(constr_poly, [], 'g-', [], confg);

hold on;

confg.plot.shade = 0.1;

plot(ellipsoid_constr, [], 'b-');

hold off;

%% function

[P_opt] = maxInscribedEllipsoid (A, b, c)

[m, n] = size(A);

% variables

P = sdpvar(n, n, 'symmetric');

% Constraints: P ≻ 0

Constraints = [P >= eps*eye(n)];

% Containment: for each i, ||P^(1/2)*a_i|| + a_i'*c <= b_i

for i = 1:m

ai = A(i, :)';

Constraints = [Constraints, norm(sqrtm(P)*ai, 2) + ai'*c <= b(i)];

end

% Objective: maximize geometric mean of eigenvalues of P

Objective = -geomean(P);

% Objective: minimize log det of P

% Objective = logdet(P);

% Solver settings: use SeDuMi or mosek

% options = sdpsettings('solver', 'sedumi', 'verbose', 1);

options = sdpsettings( 'solver', 'mosek', 'verbose', 1);

% Solve the SDP

sol = optimize(Constraints, Objective, options)

% Extract optimal values

P_opt = value(P);

end

Result:
% solvertime: 0
% info: 'No suitable solver'
% problem: -2
% yalmiptime: 0.2500

Can you help me find the problem, please?

[Johan Löfberg]

You have to perform a variable change as the current form is nastily nonlinear in P (norm of nonlinear square-root). The problem is linear in a square-root of P, and the objective is monotonic when doing that variable change, hence you should optimize over a new variable Q representing that root.

[Sunandan Adhikary]

Thanks.
Actually, I have another constraint which was linear in terms of P but would be bilinear if I consider sqrt_P as an sdpvar.

constr= [P Apm*P + Bpm*Z; ...

(Apm*P + Bpm*Z)' (1 + alpha_star_itr)*P] >= slack ; Apm and Bpm are characteristic matrices, and alpha_star_itr is the desired decay \in [-1, 0]. Z is an sdpvar as well. Would it be okay to consider sqrt_P as a variable in this case, or is there some other way?

[Johan Löfberg]

you will have to come up with a model which linearizes both the ellipsoidal constraints and the lyapunov constraint

[Sunandan Adhikary]

Yes. Trying that. Thanks for your responses and suggestions.

[Sunandan Adhikary]

I added a new constraint:

constr1 = [P sqrt_P;...

sqrt_P' eye(size(P))] ==0 (or <= since hard constraints are unsupported) and updated the above constraints using this sqrt_P in place of P with logdet objective. sdpt3 solver still returns numerical problems! Can you suggest any modification in the model?

[Johan Löfberg]

That constraint is trivially infeasible as you are saying I == 0

[Sunandan Adhikary]

I changed == to >=, and now for some systems YALMIP shows number of iterations exceeded even if I increase the default number of iterations for sdpt3.

Any corrections are needed in the model?

[Johan Löfberg]

>> test_smc
Undefined function or variable 'limits_to_cvxset'.

[Johan Löfberg]

and it smells like that model would be unbounded (your objective says make P1 large, and then you change equality to P1>=Psqrt*Psqrt which means it can let P1 tend to infinity in some direction)

On Monday, 14 July 2025 at 09:30:11 UTC+2 mesun...@gmail.com wrote:

[Sunandan Adhikary]

I used -logdet(-logdet(sqrt_inverseP1);) as the objective (ignore ~use_inverseP parts). The largest inscribable ellipsoid constraint constr_inv limits sqrt_inverseP1. Is that not enough?

Sorry for incomplete resources. Here are two sets of matrices and polytopes. For one, the solver is able to converge(system.mat); for another one(system-acc.mat), it keeps diverging, producing a bad solution.

[Johan Löfberg]

you will have to supply reproducible example. this fails/warns with many things

>> clear all
>> test_smc
Warning: Could not find appropriate function on path loading function
handle D:\workspace2\mlfcbf\limits_to_cvxset.m>limits_to_cvxset/evalSafety
> In test_smc (line 17)
domain variable existsWarning: Using 'state' to set RAND's internal state causes RAND, RANDI,
and RANDN to use legacy random number generators.  This syntax is not
recommended.  See Replace Discouraged Syntaxes of rand and randn to use
RNG to replace the old syntax.
> In rng>handleLegacyStruct (line 160)
  In rng (line 84)
  In rand_hash (line 15)
  In appendYALMIPvariables (line 35)
  In sdpvar (line 699)
  In test_smc (line 41)

alpha_star_itr =

  -0.037068607427671

Warning: Equality constraint evaluated to trivial true.
> In constraint (line 45)
  In  ==  (line 17)
  In test_smc (line 89)

[Sunandan Adhikary]

Sorry for the inconvenience.
Line 35 with  limits_to_cvxset () function is commented now. Here are the outputs when run in my system (changed == to <= in lines 89 and 93 to avoid warnings).
Need to change the mat file names in line 17 while the provided mat files are in the same path.

[Johan Löfberg]

You've got multiple issues going on

To begin with, the problem is unbounded, i.e. you are not modelling what you want (easily seen by removing the objective which leads to no issues, while keeping the objective leads to problems). 

You can for instance do

sol_safe = optimize([constr_inv,trace(sqrt_inverseP1)<=1e5], obj, ops)
log(det(value(sqrt_inverseP1)))

and no matter how large you make the artificial bound the objective becomes larger. I.e. you are relaxing something in the wrong direction (which sort of looks reasonable, you've got sqrt_inverseP2 bounded from above through the ellipsoidal stuff, but then you only have sqrt_inverseP1 bounded from below, and an objective where you make it large, thus trivially you can make it arbitrarily large)

Second, inverseP1 >= slack is not what you want most likely. You problably want inverseP1 >= slack*eye(length(invesrseP1))

Third, this should make you concerned

K>> constr_inv = [inverseP1 >= slack, sqrt_inverseP1 >= slack,...
                            sqrt_inverseP2 <= Cpm*sqrt_inverseP1*Cpm'];
Warning: Inequality constraint evaluated to trivial true. 

Doesn't make sense to add that constraint, since you already defined sqrt_inverseP2 to be  Cpm*sqrt_inverseP1*Cpm'.

[Sunandan Adhikary]

I understand. Yes, I mainly want a comparison between diagonals rather than an element-wise comparison. Thanks.
This means if I consider -logdet(sqrt_inverseP2) as the objective, then this problem should be successfully solved.
Logdet inherently presumes the matrix to be Hermitian. But I understand that a real matrix, sqrt_inverseP2, can also be used in an objective.
However, for some reason, it does not give a stable controller for system-acc. I am checking it.
Thanks again for the help.

[Johan Löfberg]

No, I don't think you want a comparison between diagonals (then you should do diag(A)>=diag(B)). You want A to be positive definite, A >= slack*eye(n). As you wrote it now, it is interpreted as A >= slack*ones(n) which is something very different

[Sunandan Adhikary]

Yes. I understood.
Now it is being successfully solved, but trying to figure out why the decay constraint does not ensure system stability!

[Sunandan Adhikary]

Positive PFEAS are basically primal residuals that ensure the satisfaction of a certain constraint, right? 

[Johan Löfberg]

primals in sdpt3 are duals in your view of the model, i.e.  feasibility of the primal  is related to optimality conditions