# How to fix such type of problem ?

3 views (last 30 days)
Assen Beshr on 29 Nov 2023
Answered: Sam Chak on 29 Nov 2023
% 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

Sam Chak on 29 Nov 2023
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