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
