LMI optimization problem

Hi there

I am regenerating an optimization algorithm from a paper and I am using exactly the example provided by the ref.

this is the algorithm:

It starts a long operation (continuously trying to obtain results - non-stop)- for example, it is +20 min trying and not stopped yet and my system (mac M2) is completely hot, I have never seen it in this way!
during several loops bring this answer:
(I am sure it must obtained and there is a solution for this example I think there is a problem with my code but I cannot find it…

``````Lin. dep.  - primal attempts        : 0                 successes              : 0
Lin. dep.  - dual attempts          : 0                 successes              : 0
Lin. dep.  - primal deps.           : 0                 dual deps.             : 0
Presolve terminated. Time: 0.00
Optimizer terminated. Time: 0.00

υ ≲ 0. A controller may not be designed using this algorithm.
Problem
Name                   :
Objective sense        : minimize
Type                   : CONIC (conic optimization problem)
Constraints            : 45
Affine conic cons.     : 0
Disjunctive cons.      : 0
Cones                  : 0
Scalar variables       : 2
Matrix variables       : 1 (scalarized: 6)
Integer variables      : 0

Optimizer started.
Presolve started.
Eliminator - tries                  : 0                 time                   : 0.00
Lin. dep.  - tries                  : 0                 time                   : 0.00
Lin. dep.  - primal attempts        : 0                 successes              : 0
Lin. dep.  - dual attempts          : 0                 successes              : 0
Lin. dep.  - primal deps.           : 0                 dual deps.             : 0
Presolve terminated. Time: 0.00
Optimizer terminated. Time: 0.00

υ ≲ 0. A controller may not be designed using this algorithm.
Problem
Name                   :
Objective sense        : minimize
Type                   : CONIC (conic optimization problem)
Constraints            : 45
Affine conic cons.     : 0
Disjunctive cons.      : 0
Cones                  : 0
Scalar variables       : 2
Matrix variables       : 1 (scalarized: 6)
Integer variables      : 0

Optimizer started.
Presolve started.
.
.
.

``````

my code: (in my code, I used epsilon for myself, and used upsilon instead of ref. epsilon, so upsilone in the code, is equal to epsilon in the screenshot

``````using JuMP
using Mosek
using Convex
using MosekTools
using LinearAlgebra
using ControlSystems
using LinearOperators
using MatrixEquations

A = [-0.15 1.90 1.55;
0.50 -0.30 0.10;
0.20 0.50 -2.55]

B = [0.55 -0.64 0.16;
1.69 0.38 0;
0.59 -1.50 1.31]

C = [1 1 0;
0 0 1]

Q = [1.841535183506899 0.6798746132078772 1.1829150444507845;
0.6798746132078772 0.5221192150230858 0.24559457200861234;
1.1829150444507845 0.24559457200861234 0.9063941504567442]

# Assisst matrices
R = B*B'

j = 1
β_j = 1
δ = 0.1

# (9): A'X + XA − XBB' X + Q = 0 ----> A'X + XA − XRX + Q = 0
X, CLSEIG = arec(A,R,2I)

# (5): [(A−BB'X)'P+P(A−BB'X),X'B; B'X,-σ1*I]<0 and (6): A'P + PA−σ2 C'C < 0
while true
n = size(A, 1)
r = size(B, 2)
ϵ = 1e-10
Xn =  β_j * X

model = Model(Mosek.Optimizer)
@variable(model, σ1 ≥ 0)
@variable(model, σ2 ≥ 0)
@variable(model, P[1:n, 1:n], PSD)

@constraint(model, [(A-R*X)'*P+P*(A-R*X) X'*B;B'*X -σ1*Matrix(I,r,r)] .<= -ϵ)
@constraint(model, A'*P+P*A-σ2*C'*C .<= -ϵ)
optimize!(model)

# Check feasibility
is_solved_and_feasible(model; allow_local = false)

# Define feasible variable
feasible = is_solved_and_feasible(model)

if feasible
# Get the optimized values of P, σ1, and σ2
P_star = value.(P)
σ1_star = value(σ1)
σ2_star = value(σ2)

return feasible, P_star, σ1_star, σ2_star
m = Model(Mosek.Optimizer)
@variable(m, K[1:n, 1:n])  # Controller gain matrix K
# Matrix inequality (10)
# Constrain decision variables (replace with actual data if available)
for i in 1:n
for j in 1:n
if i != j
@constraint(m, [A[i, j]+B[i, :]*K*C[:, j]] >= 0)
end
end
end
# Matrix inequality (11)
@constraint(m, ((A+B*K*C)'*P_star+P_star*(A+B*K*C)) <= -ϵ)
# Solve the optimization problem
optimize!(m)

# Check feasibility
is_solved_and_feasible(m; allow_local = false)

# Define feasible variable
correct = is_solved_and_feasible(m)

if correct
println("Solution found!")
println("Designed static controller K:")
else
break
end
else
υ = β_j - δ
if υ <= -ϵ
println("υ ≲ 0. A controller may not be designed using this algorithm.")
else
# Update beta_j
global  j += 1
for j in 1:10000
global β_j = β_j - δ
continue
end
end
end
end
``````

I can’t run your code because it doesn’t seem to be all there.

It has, for example, `return feasible, P_star, σ1_star, σ2_star`.

Also, after `println("υ ≲ 0. A controller may not be designed using this algorithm.")`, it appears that you don’t break, so you end up in an infinite loop.

1 Like

Dear Odow

The code is complete, this part:

``````# Check feasibility
is_solved_and_feasible(model; allow_local = false)

# Define feasible variable
feasible = is_solved_and_feasible(model)

if feasible
# Get the optimized values of P, σ1, and σ2
P_star = value.(P)
σ1_star = value(σ1)
σ2_star = value(σ2)
``````

is related to:

``````model = Model(Mosek.Optimizer)
@variable(model, σ1 ≥ 0)
@variable(model, σ2 ≥ 0)
@variable(model, P[1:n, 1:n], PSD)

@constraint(model, [(A-R*X)'*P+P*(A-R*X) X'*B;B'*X -σ1*Matrix(I,r,r)] .<= -ϵ)
@constraint(model, A'*P+P*A-σ2*C'*C .<= -ϵ)
``````

code is complete on my side, it goes to an inf loop as you said.
Should I add “break” after the mentioned part?
in that case, it cannot find “proper K”?

what is the problem when I am sure there is a solution?
I am sure the is wrong with my code.

You have a number of issues.

`@constraint(model, A'*P+P*A-σ2*C'*C .<= -ϵ)`

What do you intend this constraint to be? (There are others like it.) Currently, it constrains each element of the matrix `A'*P+P*A-σ2*C'*C` to be `<= -ϵ`.

``````@constraint(model, A'*P+P*A-σ2*C'*C <= 0, PSDCone())
``````

`A[i, j]+B[i, :]*K*C[:, j]`

This doesn’t make sense. `K` is `3 * 3`, but `C` is `2 * 3`. You cannot multiply `K * C[:, j]`.

1 Like

I replaced this part:

``````@constraint(model, A'*P+P*A-σ2*C'*C <= 0, PSDCone())
``````

and I added “break” and now, it is not inf-loop.
but bring this result:

``````"υ ≲ 0. A controller may not be designed using this algorithm."
``````

while I am sure, it must be able to find the proper K,
Because I am using the tested example, I am sure the problem is with my code.

about the size of K, we have `A + B*K*C` where A = 3 by 3, B = 3 by 3, C = 2 by 3

and we are trying to obtain proper K using the algorithm, K must be 3 by 2.

first we must do BK and results * C----> (BK)*C in that case, the sizes are true,

in MATLAB we always do this in this way, K will be found considering sizes B and C.

but here, if this is the problem, I do not how to solve it, and if this is the problem, why there is no mismatch dimension err?

• update:

I think maybe the problem is in step `V` algorithm, this section of the code:

``````    # Check feasibility
is_solved_and_feasible(m; allow_local = false)

# Define feasible variable
correct = is_solved_and_feasible(m)

if correct
println("Solution found!")
println("Designed static controller K:")
else
break
end
else
υ = β_j - δ
if υ <= -ϵ
println("υ ≲ 0. A controller may not be designed using this algorithm.")
break
else
# Update beta_j
global  j += 1
for j in 1:10000
global β_j = β_j - δ
continue
end
end
end
end
``````

if the system is controllable, and proper K could not be found in the first iteration (j=1),
it must go back to step `III`, I think I have a problem in this part (I am not sure).

we are trying to obtain proper K using the algorithm, K must be 3 by 2

Sure. Except that you have defined it as 3x3: `@variable(m, K[1:n, 1:n])`

``````global  j += 1
for j in 1:10000
global β_j = β_j - δ
continue
end
``````

What do you think this does? Why the loop 10,000 times? I think you just want instead: `global β_j = β_j - δ`.

I don’t really know if this is exactly what you want, but it should point you in the right direction:

``````julia> using JuMP

julia> using SCS

julia> using LinearAlgebra

julia> using MatrixEquations

julia> function main(A, B, C, Q)
R = B * B'
β_j = 1
δ = 0.1
n = size(A, 1)
r = size(B, 2)
X, _ = MatrixEquations.arec(A, R, 2 * LinearAlgebra.I)
while β_j > 0
Xn =  β_j * X
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, σ1 >= 0)
@variable(model, σ2 >= 0)
@variable(model, P[1:n, 1:n], PSD)
@constraint(
model,
[
((A - R * X)' * P + P * (A - R * X)) (X' * B)
(B' * X) -σ1*LinearAlgebra.I(r)
] <= 0, PSDCone())
@constraint(model, A' * P + P * A - σ2 * C' * C <= 0, PSDCone())
optimize!(model)
if !is_solved_and_feasible(model)
β_j = β_j - δ
continue
end
P_star = value.(P)
m = Model(SCS.Optimizer)
set_silent(m)
@variable(m, K[1:n, 1:size(C, 1)])
for i in 1:n, j in 1:n
if i != j
@constraint(m, A[i, j] + B[i, :]' * K * C[:, j] >= 0)
end
end
@constraint(m, (A+B*K*C)'*P_star + P_star*(A+B*K*C) <= 0, PSDCone())
optimize!(m)
if is_solved_and_feasible(m)
println("Designed static controller K:")
return value.(K)
else
break
end
end
println("υ ≲ 0. A controller may not be designed using this algorithm.")
return nothing

end
main (generic function with 1 method)

julia> A = [
-0.15 1.90 1.55
0.50 -0.30 0.10
0.20 0.50 -2.55
]
3×3 Matrix{Float64}:
-0.15   1.9   1.55
0.5   -0.3   0.1
0.2    0.5  -2.55

julia> B = [
0.55 -0.64 0.16
1.69 0.38 0
0.59 -1.50 1.31
]
3×3 Matrix{Float64}:
0.55  -0.64  0.16
1.69   0.38  0.0
0.59  -1.5   1.31

julia> C = [
1 1 0
0 0 1
]
2×3 Matrix{Int64}:
1  1  0
0  0  1

julia> Q = [
1.841535183506899 0.6798746132078772 1.1829150444507845
0.6798746132078772 0.5221192150230858 0.24559457200861234
1.1829150444507845 0.24559457200861234 0.9063941504567442
]
3×3 Matrix{Float64}:
1.84154   0.679875  1.18292
0.679875  0.522119  0.245595
1.18292   0.245595  0.906394

julia> main(A, B, C, Q)
Designed static controller K:
3×2 Matrix{Float64}:
-0.710233  -0.705731
2.99965    3.02519
3.66925    5.19916
``````
1 Like

Hi again,

I tried to learn what you used in the provided code, step by step-reading manuals related to the structures you used to introduce the variables and constraints,

I tried to check again and again and again the code and main algorithm, also the code gave us K successfully but it is not the correct controller,

Because closed-loop eigenvalues for A+BKC are:

``````2.5216
``````

-1.0481
-4.5343

after several revised, I changed 2 points:

1. step V according to the first post, I added it,
2. I replaced υ character with η because I found that this character has a conflict with sth.

so the code:

``````using JuMP
using SCS
using LinearAlgebra
using MatrixEquations

function main(A, B, C, Q)
R = B * B'
β_j = 1
δ = 0.1
n = size(A, 1)
r = size(B, 2)
X, _ = MatrixEquations.arec(A, R, Q)
while β_j > 0
Xn =  β_j * X
model = Model(SCS.Optimizer)
set_silent(model)
@variable(model, σ1 >= 0)
@variable(model, σ2 >= 0)
@variable(model, P[1:n, 1:n], PSD)
@constraint(
model,
[
((A - R * Xn)' * P + P * (A - R * Xn)) (Xn' * B)
(B' * Xn) -σ1*LinearAlgebra.I(r)
] <= 0, PSDCone())
@constraint(model, A' * P + P * A - σ2 * C' * C <= 0, PSDCone())
optimize!(model)
if !is_solved_and_feasible(model)
β_j = β_j - δ
continue
end
P_star = value.(P)
m = Model(SCS.Optimizer)
set_silent(m)
@variable(m, K[1:n, 1:size(C, 1)])
for i in 1:n, j in 1:n
if i != j
@constraint(m, A[i, j] + B[i, :]' * K * C[:, j] >= 0)
end
end
@constraint(m, (A+B*K*C)'*P_star + P_star*(A+B*K*C) <= 0, PSDCone())
optimize!(m)
if is_solved_and_feasible(m)
println("Designed static controller K:")
return value.(K)
else
η = β_j - δ
if η <= -ϵ
println("η ≲ 0. A controller may not be designed using this algorithm.")
return nothing
break
else
# Update beta_j
j = 0
while β_j > 0
j += 1
β_j = β_j - δ
end
continue
end
end
end
end
``````

but I received still the same K as the controller and it does not stabilize the sys A-B-C

in the algorithm, we have constraints <0 but to path the err, we used <=0, could it be the problem?

Would you mind sharing a reference to the paper from which the algorithm is coming?

1 Like

yes sure,

Static Output-Feedback Stabilization for MIMO LTI Positive Systems using LMI-based Iterative Algorithms

any help can save me …

1 Like

I don’t think your beta update matches your screenshot. It should be updated only once each iteration. Your code currently updates it until it is negative.