Problem with minimizing a LQR problem with Genetic Algorithm

I am trying to solve a LQR problem using Genetic Algorithm.

The values of the LQR matrices are not importent to me, I am just trying to understand the way that it is work.

First I solve the system using the Matlab's LQR algorithm and them I try to solve using the Genetic Algorithm and compare the results.

In the code attach I am trying to compare between P to P2, but I dont get a match.

Also when I run the code, I get the following Error: Optimization terminated: average change in the fitness value less than options

%% State Space
A = [1 2; 3 4];
B = [1; 1];
Q = [1 0;0 1];
R = 1;

%% Slove the LQR
[K,P1] = lqr(A,B,Q,R);
P2 = norm(P1);

%% Solve the Genetic Algorithm
P = ga(@cost_fun,2);

Optimization terminated: maximum number of generations exceeded.

function P = cost_fun(K1)
A = [1 2; 3 4];
B = [1; 0];
Q = [1 0;0 1];
R = 1;
A_i = A - B*K1;
Q_i = - Q - K1'*R*K1;
P_i = lyap(A_i,Q_i);
P = norm(P_i);
end

Answer 1

[Sam Chak]

Hi @gilad shaul

Minimizing the norm(P) does not work because there are infinite K that can make A'*P + P*A - K'*R*K + Q exactly 0. If you want to solve using genetic algorithm, ga(), then you should directly minimize the Finite-time Performance Index J (a scalar) instead:

Else, if you want to find the elements in the positive-definite matrix that minimizes J, then you should use gamultiobj() instead (because there are multiple fitness functions). In 2nd-order systems, there are three fitness functions involving , , because the positive-definite matrix looks like this:

.

% ------- Computation of K via LQR (solving Algebraic Riccati Matrix Equation) -------
% State Space
A = [0 1; 0 -1];
B = [0; 1];
Q = [1 0; 0 1];
R = 1;

% LQR
K = lqr(A, B, Q, R)

K = 1×2

1.0000 1.0000

% [K, P, E] = lqr(A, B, Q, R)
% Riccati Equation
% A'*P + P*A - (P*B)*inv(R)*(B'*P) + Q
% QQ = - (P*B)*inv(R)*(B'*P) + Q;
% QQ = - (P*B)*K + Q;
% QQ = - K'*R*K + Q;
% PP = lyap(A', QQ)     % returns the same positive-definite matrix P

% ------- Computation of K via GA (minimizing Finite-time Performance Index J) -------
fitfun  = @costfun;
PopSize = 100;
MaxGens = 200;
nvars   = 2;
A       = -eye(nvars);
b       = zeros(nvars, 1);
Aeq     = [];               
beq     = [];
lb      = [];
ub      = [];
nonlcon = [];

options = optimoptions('ga', 'PopulationSize', PopSize, 'MaxGenerations', MaxGens);
[KK, fval, exitflag, output] = ga(fitfun, nvars, A, b, Aeq, beq, lb, ub, nonlcon, options)

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

KK = 1×2

0.9986 1.0004

fval = 2.0024

exitflag = 1

output = struct with fields:
      problemtype: 'linearconstraints'
         rngstate: [1×1 struct]
      generations: 56
        funccount: 5425
          message: 'Optimization terminated: average change in the fitness value less than options.FunctionTolerance.'
    maxconstraint: 0
       hybridflag: []

function J = costfun(param)
    % Parameters
    A      = [0 1; 0 -1];
    B      = [0; 1];
    Q      = [1 0; 0 1];
    R      = 1;
    K      = [param(1) param(2)];
    
    % State-Space
    odefcn = @(t, x) (A - B*K)*x;
    tspan  = [0 20];    % Finite-time
    x0     = [1 0];
    [t, x] = ode45(odefcn, tspan, x0);
    
    % Performance Index
    u      = - K(1)*x(:,1) - K(2)*x(:,2);
    itgrd  = Q(1,1)*x(:,1).^2 + Q(2,2)*x(:,2).^2 + R*u.^2;  % integrand
    J = trapz(t, itgrd);
end

[Gilad Shaul]

Thank you very much!

That was really helpful.

If I want to find the K that minimize P, how do you suggest to do that?

Answer 2

[Sam Chak]

Hi @gilad shaul

To find through , you need to expand the Continuous-time Algebraic Riccati Equation (CARE) in terms of , and then use gamultiobj() to solve the root-finding problem.

If you like the alternative solution, please consider voting 👍 the Answer as a little token of appreciation.

%% ----- Expanding the CARE -----
syms p11 p12 p22

A = sym('A', [2 2]);	% state matrix
B = sym('B', [2 1]);	% input matrix
P = sym('P', [2 2]);	% positive definite matrix

A = [sym('0') sym('1');
     sym('0') sym('-1')];

B = [sym('0'); sym('1')];

P = [p11 p12;
     p12 p22];

Q = sym(eye(2));

R = sym('1');

CARE = A.'*P + P*A - P*B/R*B.'*P + Q

CARE = 

%% ----- Solving using gamultiobj -----
fitfun  = @multiobjfcn;
PopSize = 15;               
MaxGen  = 30;
nvars   = 3;                
A       = -eye(nvars);      
b       = zeros(nvars, 1);  
Aeq     = [];               
beq     = [];               
lb      = [0 0 0];          
ub      = [3 3 3];          
nonlcon = [];               
options = optimoptions('gamultiobj', 'ConstraintTolerance', 1e-6, 'FunctionTolerance', 1e-6, 'HybridFcn', @fgoalattain);

[p, fval] = gamultiobj(fitfun, nvars, A, b, Aeq, beq, lb, ub, nonlcon, options); 

Optimization terminated: maximum number of generations exceeded.

popt = p(end, :);
R = 1;
B = [0; 1];
P = [popt(1) popt(2); popt(2) popt(3)]

P = 2×2

2.0000 1.0000 1.0000 1.0000

Kopt = R\(B'*P)

Kopt = 1×2

1.0000 1.0000

%% ----- Multi-objective Functions -----
function y = multiobjfcn(p)
    y    = zeros(3, 1);
    y(1) = (1 - p(2)^2)^2;                      % Remember to square the CARE
    y(2) = (- p(3)^2 - 2*p(3) + 2*p(2) + 1)^2;
    y(3) = (p(1) - p(2) - p(2)*p(3))^2;
end