基于遗传算法的LQR控制器最优设计算法

线性二次调节器(LQR)是控制理论中重要的设计方法，而遗传算法(GA)为LQR控制器的优化设计提供了强大的全局搜索能力。

LQR控制器基本原理

LQR控制器通过最小化代价函数设计最优状态反馈增益矩阵：

J=∫(xTQx+uTRu)dtJ = ∫(xᵀQx + uᵀRu)dtJ=∫(xTQx+uTRu)dt

其中：

• Q≥0Q ≥ 0Q≥0：状态加权矩阵
• R>0R > 0R>0：控制输入加权矩阵
• u=−Kxu = -Kxu=−Kx：状态反馈控制律

传统方法通过求解代数RiccatiRiccatiRiccati方程获得最优增益矩阵KKK：

K=R−1BTPK = R⁻¹BᵀPK=R−1BTP
AP+PAT−PBR−1BTP+Q=0AP + PAᵀ - PBR⁻¹BᵀP + Q = 0AP+PAT−PBR−1BTP+Q=0

遗传算法优化LQR设计

优化问题分析

1. 优化变量：Q和R矩阵的元素
2. 目标函数：闭环系统性能指标（如响应时间、超调量、控制能耗）
3. 约束条件：系统稳定性、控制输入限制

算法流程

否

是

初始化种群

评估适应度

满足终止条件?

选择操作

交叉操作

变异操作

新一代种群

输出最优解

MATLAB实现

1. 系统建模与适应度函数

function[fitness, stability]=lqr_fitness_function(params, A, B, Q_structure, R_structure)% 提取参数 n =size(A,1); m =size(B,2);% 重构Q矩阵ifstrcmp(Q_structure,'diagonal') q_diag =params(1:n); Q =diag(q_diag); param_offset = n;elseifstrcmp(Q_structure,'full') Q =reshape(params(1:n*n), n, n); Q =(Q + Q')/2;% 确保对称性 param_offset = n*n;end% 重构R矩阵ifstrcmp(R_structure,'scalar') R =params(param_offset +1); R = R *eye(m);elseifstrcmp(R_structure,'diagonal') r_diag =params(param_offset +1:param_offset + m); R =diag(r_diag); param_offset = param_offset + m;elseifstrcmp(R_structure,'full') R =reshape(params(param_offset +1:param_offset + m*m), m, m); R =(R + R')/2;% 确保对称性end% 求解Riccati方程try[K,~, eig_vals]=lqr(A, B, Q, R); stability =all(real(eig_vals)<0);% 检查稳定性catch fitness =inf;% 求解失败，赋予极大代价 stability = false;return;end% 闭环系统仿真 sys_cl =ss(A - B*K, B,eye(n),0); t =0:0.01:10; x0 =[1;zeros(n-1,1)];% 初始状态[y,~, x]=initial(sys_cl, x0, t);% 计算性能指标 settling_time =settling_time_calc(y(:,1), t);% 调节时间 overshoot =max(y(:,1))/y(1,1)-1;% 超调量 control_effort =sum(sum((K*x).^2))*0.01;% 控制能耗% 加权适应度函数 w1 =0.4;% 调节时间权重 w2 =0.3;% 超调量权重 w3 =0.3;% 控制能耗权重 fitness = w1*(settling_time/10)+ w2*overshoot + w3*(control_effort/100);% 惩罚不稳定系统if~stability fitness = fitness +10;endendfunction ts =settling_time_calc(response, time) final_value =response(end); tolerance =0.02*abs(final_value);% 2%容差带 idx =find(abs(response - final_value)<= tolerance,1,'first');ifisempty(idx) ts =time(end);else ts =time(idx);endend

2. 遗传算法主程序

function[best_K, best_params, best_fitness]=ga_lqr_design(A, B, Q_structure, R_structure)% 参数设置 n =size(A,1);% 状态维度 m =size(B,2);% 输入维度% 根据矩阵结构确定变量个数ifstrcmp(Q_structure,'diagonal') nQ = n;elseifstrcmp(Q_structure,'full') nQ = n*n;endifstrcmp(R_structure,'scalar') nR =1;elseifstrcmp(R_structure,'diagonal') nR = m;elseifstrcmp(R_structure,'full') nR = m*m;end total_vars = nQ + nR;% 遗传算法参数 pop_size =50; max_gen =100; crossover_rate =0.8; mutation_rate =0.1;% 变量边界 (Q对角元素>0, R元素>0) lb =zeros(total_vars,1); ub =ones(total_vars,1)*100;% 初始化种群 population = lb +(ub - lb).*rand(pop_size, total_vars);% 进化循环 best_fitness_history =zeros(max_gen,1);for gen =1:max_gen % 评估适应度 fitness =zeros(pop_size,1); stability_flags =false(pop_size,1);parfori=1:pop_size [fit, stable]=lqr_fitness_function(population(i,:), A, B, Q_structure, R_structure);fitness(i)= fit;stability_flags(i)= stable;end% 记录最佳个体[best_fit, idx]=min(fitness); best_individual =population(idx,:);best_fitness_history(gen)= best_fit;fprintf('Generation %d: Best Fitness = %.4f\n', gen, best_fit);% 选择操作 (锦标赛选择) new_population =zeros(size(population));fori=1:pop_size candidates =randperm(pop_size,2);[~, winner]=min(fitness(candidates));new_population(i,:)=population(candidates(winner),:);end% 交叉操作 (算术交叉)fori=1:2:pop_size-1if rand < crossover_rate alpha = rand; parent1 =new_population(i,:); parent2 =new_population(i+1,:); child1 = alpha*parent1 +(1-alpha)*parent2; child2 = alpha*parent2 +(1-alpha)*parent1;new_population(i,:)= child1;new_population(i+1,:)= child2;endend% 变异操作 (高斯变异)fori=1:pop_size if rand < mutation_rate mutation =randn(1, total_vars)*0.1;new_population(i,:)=new_population(i,:)+ mutation;% 确保在边界内new_population(i,:)=max(min(new_population(i,:), ub), lb);endend population = new_population;end% 提取最优解 best_params = best_individual;[~, stability]=lqr_fitness_function(best_params, A, B, Q_structure, R_structure);% 重构最优Q和Rifstrcmp(Q_structure,'diagonal') q_diag =best_params(1:n); Q_opt =diag(q_diag); param_offset = n;elseifstrcmp(Q_structure,'full') Q_opt =reshape(best_params(1:n*n), n, n); Q_opt =(Q_opt + Q_opt')/2; param_offset = n*n;endifstrcmp(R_structure,'scalar') R_opt =best_params(param_offset +1)*eye(m);elseifstrcmp(R_structure,'diagonal') r_diag =best_params(param_offset +1:param_offset + m); R_opt =diag(r_diag);elseifstrcmp(R_structure,'full') R_opt =reshape(best_params(param_offset +1:end), m, m); R_opt =(R_opt + R_opt')/2;end% 计算最优增益矩阵 best_K =lqr(A, B, Q_opt, R_opt); best_fitness = best_fit;% 绘制进化过程 figure;plot(1:max_gen, best_fitness_history,'LineWidth',2);xlabel('Generation');ylabel('Best Fitness');title('Genetic Algorithm Convergence'); grid on;end

3. 示例使用案例

% 示例系统：倒立摆 m =0.5;% 小车质量 (kg) M =0.5;% 摆杆质量 (kg) l =0.3;% 摆杆长度 (m) g =9.81;% 重力加速度 (m/s²)% 状态空间矩阵 A =[0100;00-(m*g)/M 0;0001;00((M+m)*g)/(M*l)0]; B =[0;1/M;0;-1/(M*l)];% 使用遗传算法优化LQR参数% Q结构：对角矩阵 (状态加权)% R结构：标量 (控制输入加权)[Q_struct, R_struct]=deal('diagonal','scalar');[best_K, best_params, best_fitness]=ga_lqr_design(A, B, Q_struct, R_struct);% 显示结果disp('Optimal gain matrix K:');disp(best_K);% 重构最优Q和R n =size(A,1);ifstrcmp(Q_struct,'diagonal') q_diag =best_params(1:n); Q_opt =diag(q_diag); R_opt =best_params(n+1);elseifstrcmp(Q_struct,'full') Q_opt =reshape(best_params(1:n*n), n, n); Q_opt =(Q_opt + Q_opt')/2; R_opt =best_params(n*n+1);end% 验证性能 sys_open =ss(A, B,eye(4),0); sys_closed =ss(A - B*best_K, B,eye(4),0);% 仿真对比 t =0:0.01:10; x0 =[0.1;0;0.1;0];% 初始状态 figure;subplot(2,1,1);initial(sys_open, x0, t);title('Open-loop Response');legend('Position','Velocity','Angle','Angular Velocity');subplot(2,1,2);initial(sys_closed, x0, t);title('Closed-loop Response with GA-LQR');legend('Position','Velocity','Angle','Angular Velocity');% 性能指标对比[y_open, t_open]=initial(sys_open, x0, t);[y_closed, t_closed]=initial(sys_closed, x0, t); overshoot_open =max(y_open(:,3))-y_open(1,3); overshoot_closed =max(y_closed(:,3))-y_closed(1,3);fprintf('Open-loop overshoot: %.4f rad\n', overshoot_open);fprintf('GA-LQR closed-loop overshoot: %.4f rad\n', overshoot_closed);

算法改进策略

1. 混合优化策略

function[K_opt, params_opt]=hybrid_optimization(A, B)% 第一阶段：遗传算法粗搜索[~, params_ga,~]=ga_lqr_design(A, B,'diagonal','scalar');% 第二阶段：局部搜索精调 options =optimoptions('fmincon','Algorithm','sqp',...'Display','iter','StepTolerance',1e-6); fun =@(x)lqr_fitness_function(x, A, B,'diagonal','scalar'); params_opt =fmincon(fun, params_ga,[],[],[],[],...zeros(size(params_ga)),inf(size(params_ga)),...[], options);% 重构最优Q和R n =size(A,1); q_diag =params_opt(1:n); Q_opt =diag(q_diag); R_opt =params_opt(n+1);% 计算最优增益 K_opt =lqr(A, B, Q_opt, R_opt);end

2. 多目标优化

function multi_objective_fitness =multi_obj_fitness(params, A, B)% 重构Q和R% ... (同上)% 求解Riccati方程[K,~, eig_vals]=lqr(A, B, Q, R); stability =all(real(eig_vals)<0);% 闭环系统仿真 sys_cl =ss(A - B*K, B,eye(size(A,1)),0); t =0:0.01:10; x0 =rand(size(A,1),1);% 随机初始状态[y,~, x]=initial(sys_cl, x0, t);% 目标1：调节时间 settling_time =settling_time_calc(y(:,1), t);% 目标2：控制能耗 control_energy =sum(sum((K*x).^2))*0.01;% 目标3：鲁棒性 (H∞范数近似) robustness =norm(ss(A-B*K, B,eye(size(A,1)),0),'inf');% 多目标适应度 (归一化处理) multi_objective_fitness =[settling_time/max_time,... control_energy/max_energy,... robustness/max_robustness];if~stability multi_objective_fitness = multi_objective_fitness +10;endend

性能评估与可视化

functionplot_performance_comparison(A, B, K_ga, K_classic)% 测试不同初始条件 test_cases ={[0.1;0;0.1;0],[0.2;0;0.2;0],[0.3;0;0.3;0]}; t =0:0.01:5; figure;fori=1:length(test_cases) x0 = test_cases{i};% 经典LQR (Q=I, R=1) K_c =lqr(A, B,eye(4),1); sys_c =ss(A - B*K_c, B,eye(4),0);[y_c,~]=initial(sys_c, x0, t);% GA优化LQR sys_ga =ss(A - B*K_ga, B,eye(4),0);[y_ga,~]=initial(sys_ga, x0, t);% 绘制角度响应subplot(2,2,i);plot(t,y_c(:,3),'b--', t,y_ga(:,3),'r-','LineWidth',1.5);title(sprintf('Initial Angle: %.1f°',rad2deg(x0(3))));xlabel('Time (s)');ylabel('Pendulum Angle (rad)');legend('Classic LQR','GA-LQR','Location','best'); grid on;end% 性能指标对比 metrics ={'Settling Time','Overshoot','Control Effort'}; classic_metrics =zeros(3,length(test_cases)); ga_metrics =zeros(3,length(test_cases));fori=1:length(test_cases) x0 = test_cases{i};% Classic LQR sys_c =ss(A - B*K_c, B,eye(4),0);[y_c, t_c]=initial(sys_c, x0, t);classic_metrics(1,i)=settling_time_calc(y_c(:,3), t_c);classic_metrics(2,i)=max(y_c(:,3))/y_c(1,3)-1;classic_metrics(3,i)=sum(sum((K_c*initial(sys_c, x0, t_c)).^2))*0.01;% GA-LQR sys_ga =ss(A - B*K_ga, B,eye(4),0);[y_ga, t_ga]=initial(sys_ga, x0, t);ga_metrics(1,i)=settling_time_calc(y_ga(:,3), t_ga);ga_metrics(2,i)=max(y_ga(:,3))/y_ga(1,3)-1;ga_metrics(3,i)=sum(sum((K_ga*initial(sys_ga, x0, t)).^2))*0.01;end% 显示表格 figure;subplot(1,3,1);bar([mean(classic_metrics(1,:));mean(ga_metrics(1,:))]);set(gca,'XTickLabel',{'Classic LQR','GA-LQR'});title(metrics{1});ylabel('Time (s)');subplot(1,3,2);bar([mean(classic_metrics(2,:));mean(ga_metrics(2,:))]*100);set(gca,'XTickLabel',{'Classic LQR','GA-LQR'});title(metrics{2});ylabel('Percentage (%)');subplot(1,3,3);bar([mean(classic_metrics(3,:));mean(ga_metrics(3,:))]);set(gca,'XTickLabel',{'Classic LQR','GA-LQR'});title(metrics{3});ylabel('Energy Units');end

参考代码 基于遗传算法的LQR控制器最优设计算法 www.youwenfan.com/contentcsr/100888.html

应用场景与优势

典型应用领域

1. 机器人控制：机械臂轨迹跟踪、平衡控制
2. 航空航天：飞行器姿态控制、自动驾驶仪
3. 电力系统：发电机励磁控制、HVDC输电控制
4. 汽车工程：主动悬架系统、电子稳定控制

算法优势

实际应用建议

1. 问题简化：先尝试对角Q/R矩阵降低维度
2. 并行计算：利用MATLAB Parallel Computing Toolbox加速评估
3. 混合策略：结合梯度下降法进行局部搜索
4. 实时调整：设计在线自适应机制应对参数变化