How to fix such type of problem ?

% Define the objective function for LQR cost
function cost = lqr_cost(QR, A, B)
    Q = QR(1:2, :);
    R = QR(3, :);
    K = lqr(A, B, Q, R);
    eigvals = eig(A - B * K);
    cost = max(real(eigvals))^2;  % Maximize the real part of eigenvalues
end



% Define the parameters for PSO
num_particles = 20;
num_iterations = 100;
lb = [0.1*ones(2,2), 0.1*ones(1,2)]; % Lower bounds for Q and R
ub = [10*ones(2,2), 10*ones(1,2)];   % Upper bounds for Q and R

% Initialize particles with random Q and R values
particles = repmat(lb, num_particles, 1) + rand(num_particles, 5) .* (ub - lb);
velocities = zeros(num_particles, 5);
pbest = particles;
pbest_cost = inf(num_particles, 1);
gbest = zeros(1, 5);
gbest_cost = inf;

% Linear system matrices (modify A and B according to your system)
A = [0 1; -1 -1];
B = [0; 1];

% PSO optimization loop
for iter = 1:num_iterations
    for i = 1:num_particles
        % Evaluate cost for each particle
        cost = lqr_cost(reshape(particles(i, :), 2, 3), A, B);
        
        % Update personal best
        if cost < pbest_cost(i)
            pbest_cost(i) = cost;
            pbest(i, :) = particles(i, :);
        end
        
        % Update global best
        if cost < gbest_cost
            gbest_cost = cost;
            gbest = particles(i, :);
        end
    end
    
    % Update particle velocities and positions
    w = 0.7; % Inertia weight
    c1 = 1.5; % Cognitive parameter
    c2 = 1.5; % Social parameter
    r1 = rand(num_particles, 5);
    r2 = rand(num_particles, 5);
    
    velocities = w * velocities + c1 * r1 .* (pbest - particles) + c2 * r2 .* (gbest - particles);
    particles = particles + velocities;
    
    % Ensure particles stay within bounds
    particles = max(particles, lb);
    particles = min(particles, ub);
end

% Display the optimized Q and R matrices
disp('Optimized Q and R matrices:');
disp(reshape(gbest, 2, 3));

Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in PSO (line 4)

lb = [0.1*ones(2,2), 0.1*ones(1,2)]; % Lower bounds for Q and R

Answer 1

[Sam Chak]

Hi @Assen Beshr

Correct me if I interpreted your problem incorrectly. If you want to maximize the real part of the stabilizing eigenvalues (heading towards ) determined from the LQR algorithm, which requires finding the values of Q and R weights in the range , then I arrive at this result using particleswarm() optimizer:

objfun    = @costfun;
nvars     = 3;
lb        = [0.1 0.1 0.1];
ub        = [10. 10. 10.];
nonlcon   = [];
[K, fval] = particleswarm(objfun, nvars, lb, ub)

Optimization ended: relative change in the objective value 
over the last OPTIONS.MaxStallIterations iterations is less than OPTIONS.FunctionTolerance.

K = 1×3

0.1000 10.0000 0.1000

fval = 0.1421

%% Check eigenvalues
A       = [0 1; -1 -1];
B       = [0; 1];
Q       = [K(1) 0; 0 K(2)];
R       = K(3);           
K       = lqr(A, B, Q, R);
eigvals = eig(A - B*K)      

eigvals = 2×1

-0.1421 -9.9489

%% Cost function
function J = costfun(param)
    A       = [0 1; -1 -1];
    B       = [0; 1];
    Q       = [param(1) 0; 0 param(2)];
    R       = param(3);
    K       = lqr(A, B, Q, R);
    eigvals = eig(A - B*K);
    realeig = sortrows(real(eigvals));
    J       = - realeig(2);     % Maximize the real part of eigenvalues
end