# LMI optimization problem

**URL:** <https://discourse.julialang.org/t/lmi-optimization-problem/112218>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [March 28, 2024, 8:28am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218 "2024-03-28T08:28:40Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sharara](https://avatars.discourse-cdn.com/v4/letter/s/6bbea6/32.png) [@Sharara](https://discourse.julialang.org/u/Sharara)\
**Post date:** [March 28, 2024, 8:28am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/1 "2024-03-28T08:28:40Z")

</div>

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:

 ![Screenshot 0006-03-28 at 17.21.37](https://global.discourse-cdn.com/julialang/original/3X/e/d/eda54bb33b6d0ceab9b14ea9e64f70257e7ae4f9.png)

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…

```julia
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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 28, 2024, 10:54pm UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/2 "2024-03-28T22:54:14Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Sharara](https://avatars.discourse-cdn.com/v4/letter/s/6bbea6/32.png) [@Sharara](https://discourse.julialang.org/u/Sharara)\
**Post date:** [March 29, 2024, 2:23am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/3 "2024-03-29T02:23:15Z")

</div>

Dear Odow

The code is complete, this part:

```julia
# 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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 29, 2024, 4:51am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/4 "2024-03-29T04:51:08Z")

</div>

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 `<= -ϵ`.

Do you mean instead:

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

```

See [Constraints · JuMP](https://jump.dev/JuMP.jl/stable/manual/constraints/#Semidefinite-constraints)

> `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]`.

---

<div class="post-metadata">

**Author:** ![Sharara](https://avatars.discourse-cdn.com/v4/letter/s/6bbea6/32.png) [@Sharara](https://discourse.julialang.org/u/Sharara)\
**Post date:** [March 29, 2024, 11:27am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/5 "2024-03-29T11:27:00Z")

</div>

Hi Odow and Tnx for your reply.

I replaced this part:

```julia
@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:

```julia
"υ ≲ 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 B_K and results \* C----\> (B_K)\*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?

 ![Learning Julia](https://global.discourse-cdn.com/julialang/original/3X/a/5/a598545cff589e819fbb96d4a1db613473801fda.jpeg)

- update:

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

```julia
    # 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).

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 29, 2024, 9:18pm UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/6 "2024-03-29T21:18:49Z")

</div>

> 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])`

> ```julia
> 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
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

```

---

<div class="post-metadata">

**Author:** ![Sharara](https://avatars.discourse-cdn.com/v4/letter/s/6bbea6/32.png) [@Sharara](https://discourse.julialang.org/u/Sharara)\
**Post date:** [May 7, 2024, 9:30am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/7 "2024-05-07T09:30:32Z")

</div>

> [@odow](#):
>
> `main(A, B, C, Q)`

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+B_K_C 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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [May 7, 2024, 11:23am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/8 "2024-05-07T11:23:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Sharara](https://avatars.discourse-cdn.com/v4/letter/s/6bbea6/32.png) [@Sharara](https://discourse.julialang.org/u/Sharara)\
**Post date:** [May 7, 2024, 11:35am UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/9 "2024-05-07T11:35:06Z")

</div>

> [@odow](#):
>
> ` main(A, B, C, Q)`

yes sure,

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

any help can save me …

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [May 7, 2024, 9:55pm UTC](https://discourse.julialang.org/t/lmi-optimization-problem/112218/10 "2024-05-07T21:55:26Z")

</div>

> [@Sharara](#):
>
> ```julia
> j = 0
> while β_j > 0
> j += 1
> β_j = β_j - δ
> end
> 
> ```

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.