简介：这是一份基于Matlab遗传算法的LQR控制器优化设计源码包，主要面向计算机、电子信息工程及数学等专业的学生，适合用于课程设计、期末大作业或毕业设计中的控制算法仿真参考。压缩包共3个文件，包含2个m脚本和1个mdl模型文件，m脚本用于实现遗传算法优化与LQR控制器主程序，mdl文件则提供Simulink仿真模型，整体大小为93KB，结构简洁便于快速运行与二次修改。目前已有710人学习下载。借助这份源码，读者可直观了解如何通过遗传算法对LQR权重矩阵进行寻优，并搭建主动悬架等控制系统的仿真环境；代码注释清晰、可读性较好，适合有一定Matlab和自动控制基础的学习者作为参考，在理解算法流程后自行扩展功能或调整参数。
1. 为什么调LQR权重矩阵比解Riccati方程更麻烦
一个倒立摆模型摆在面前，状态方程写出来，K = lqr(A,B,Q,R)一行代码就能算出反馈增益，仿真曲线也稳定，但真正用在被控对象上却总差一口气：超调偏大、控制量饱和、某个状态收敛太慢。问题基本都出在Q和R这两个权重矩阵上。LQR 的代价函数里，Q惩罚状态偏差，R惩罚控制能量，它们的大小比例直接决定控制器是"敢动"还是"保守"，可这两个矩阵本身没有直观的物理标定方法，手调一次要反复看响应曲线，一个四阶系统调上一整天也正常。遗传算法的思路是把权重矩阵编码成个体，用种群搜索代替人工试凑，适应度函数直接取时域指标，比如 ITAE、超调量、控制能量，跑完一轮就能拿到一组接近最优的Q、R。这类基于 Matlab 的遗传算法 LQR 优化设计工程，核心代码量不大，通常就是主程序加适应度函数两段，适合做控制器参数整定、毕设预研或者工程上的快速选型。
2. 从LQR到优化目标：把控制需求写成适应度函数
2.1 LQR控制器设计的基本输入和常见痛点
LQR 设计的前提是线性时不变系统状态空间模型：
[ \dot{x} = A x + B u ]
控制目标是让二次型性能指标最小：
[ J = \int_0^\infty \left( x^T Q x + u^T R u \right) dt ]
其中Q是半正定对称矩阵，R是正定对称矩阵。Matlab 里直接用lqr(A,B,Q,R)就能求解代数黎卡提方程，得到状态反馈增益K。这里的代数计算没有任何歧义，给一组Q、R就有一组确定的K。麻烦的是Q、R的取值空间太大：Q矩阵的自由度是 (n(n+1)/2)，四阶系统就是 10 个变量，加上R，手调根本不可能全局搜索。
2.1.1 为什么Q和R不能随便给
很多人会犯的一个错误是认为Q取单位阵、R取 1 就能跑，实际结果往往是"数值上稳定，性能上平庸"。Q中对角元素代表每个状态变量的相对重要程度，比如倒立摆里角度误差比位移误差更敏感，那角度对应的Q元素就要大；R元素则代表控制输入的能量代价，R太小会让u接近执行机构极限，太大则响应迟钝。更隐蔽的是，不同状态变量的量纲不同，角度是弧度，位移是米，直接拿原始数值比较完全没意义。所以优化Q、R的第一步是理解它们在闭环极点移动中的角色，第二步才是用算法搜索。
2.2 用遗传算法搜索Q、R权重矩阵的可行编码
遗传算法要求的输入是一个实数向量，而 LQR 要求的是两个矩阵，因此要做编码变换。最常用的是把Q限制为对角矩阵，R也限制为对角矩阵，个体向量写成：
[q1, q2, ..., qn, r1, r2, ..., rm]这样四阶系统单输入时只有 5 个变量，搜索难度大幅降低。对角矩阵的合理性在于：多数控制需求只关心每个状态本身的惩罚，状态之间的交叉耦合可以在设计上放到前馈或观测器里处理，直接搜索全矩阵的自由度会让种群规模爆炸。
2.2.1 变量缩放与边界约束
遗传算法在实数编码下对变量尺度敏感。Q的元素可能从 (10^{-3}) 到 (10^3) 跨越几个量级，如果直接线性搜索，算法会花大量代数在无效区域。常见做法是对个体做对数变换，搜索空间里存的是 (\lg q_i)，解码时再用10^x还原。这样既有边界约束，又能覆盖宽量程。
变量下界和上界的设置可以参考开环特性的自然频率。下面是一组常用的边界设置：
| 变量 | 含义 | 下界 | 上界 | 说明 |
|---|---|---|---|---|
| lg(q1) | 状态1权重对数 | -4 | 4 | 对应权重 1e-4 到 1e4 |
| lg(q2) | 状态2权重对数 | -4 | 4 | 同上 |
| lg(r1) | 控制权重对数 | -6 | 2 | 控制代价通常比状态权重小 |
边界不一定要对称，如果控制量特别珍贵，R的搜索范围可以整体上移。初始种群使用PopInitRange指定，范围不要超出边界。
2.3 适应度函数怎么设才不白跑
适应度函数决定了遗传算法在朝哪个方向进化。LQR 的代价函数本身就是一种适应度，但它不包含时域约束，优化结果可能出现理论最优但超调 40% 的情况。更实用的做法是在 LQR 代价基础上叠加时域指标，比如 ITAE（时间乘绝对误差积分）或 IAE，再乘一个控制能量惩罚项。
function score = lqr_fitness(vars, sys, x0, dt, tfinal) n = size(sys.A, 1); q_vec = 10.^vars(1:n); r_vec = 10.^vars(n+1:end); Q = diag(q_vec); R = diag(r_vec); try K = lqr(sys.A, sys.B, Q, R); catch score = 1e8; return; end Acl = sys.A - sys.B * K; x = x0; t = 0:dt:tfinal; integral = 0; u_norm = 0; for k = 2:length(t) u = -K * x; x = x + dt * (Acl * x + sys.B * u); integral = integral + dt * abs(x(1)) * t(k); u_norm = u_norm + dt * u^2; end score = integral + 0.01 * u_norm; end这段代码的关键点有三个。第一，用try-catch包住lqr调用，一旦Q、R不满足正定条件或者系统不可镇定，lqr会报错，此时直接给一个超大惩罚值，让该个体被淘汰。第二，没有真的用lsim，而是用一阶欧拉递推做闭环仿真，因为遗传算法要跑几百代，每代几十个个体，每次都调用lsim会非常慢，欧拉递推在控制周期足够小时精度足够。第三，integral里乘上了t(k)，这就是 ITAE，它对末期误差不敏感，更适合工程上"早点稳定"的需求。u_norm前面的系数 0.01 是惩罚权重，如果控制量总在饱和，就调大这个系数。
3. 用Matlab实现遗传算法与LQR联合仿真
3.1 主程序框架：GA工具箱与状态方程的接口
Matlab 自带的全局优化工具箱提供了ga和gamultiobj函数，不需要额外安装第三方包。主程序要做的就是把被控对象的状态矩阵、初始条件、仿真参数准备好，然后调用ga，最后把最优个体解码成Q、R和K。
% ga_lqr_main.m rng(42); A = [0 1 0 0; 0 0 -9.8 0; 0 0 0 1; 0 0 19.6 0]; % 倒立摆线性化模型 B = [0; 1; 0; 1]; sys = ss(A, B, eye(4), 0); nvars = 5; % q1,q2,q3,q4,r1 lb = [-4, -4, -4, -4, -6]; ub = [4, 4, 4, 4, 2]; PopInitRange = [lb; ub]; options = optimoptions('ga', ... 'PopulationSize', 60, ... 'MaxGenerations', 80, ... 'CrossoverFraction', 0.8, ... 'PopInitRange', PopInitRange, ... 'Display', 'iter', ... 'UseParallel', true); fitness_fcn = @(vars) lqr_fitness(vars, sys, [0.1;0;0.05;0], 1e-3, 5); [vars_opt, fval] = ga(fitness_fcn, nvars, [], [], [], [], lb, ub, [], options); n = size(A, 1); Q_opt = diag(10.^vars_opt(1:n)); R_opt = diag(10.^vars_opt(n+1:end)); K_opt = lqr(A, B, Q_opt, R_opt); disp('最优Q对角元素：'); disp(diag(Q_opt)); disp('最优R对角元素：'); disp(diag(R_opt)); disp('LQR增益K：'); disp(K_opt);这段主程序里，rng(42)固定随机种子，保证每次运行结果一致，否则遗传算法的初始种群是随机的，同样参数可能跑出不同结果，不利于对比实验。nvars=5对应四个状态权重和一个控制权重。UseParallel设为true可以在多核 CPU 上并行计算适应度，但需要提前用parpool开启并行池，否则这个选项会被忽略。这里被控对象的A、B是倒立摆在平衡点的线性化结果，换成自己的模型时，只需要改A、B和初值。
3.2 关键参数配置：种群、交叉变异、边界约束
遗传算法的效果很大程度卡在参数选择上。ga的默认参数对控制器优化这类中等规模问题并不理想，建议按下面的表调整：
| 参数名 | 建议值 | 作用 |
|---|---|---|
| PopulationSize | 40~80 | 太小容易早熟，太大收敛慢 |
| MaxGenerations | 60~120 | 每一代耗时和收敛判断的平衡 |
| CrossoverFraction | 0.7~0.9 | 交叉比例高，种群多样性好 |
| MutationFcn | @mutationadaptfeasible | 自适应变异，适合带边界约束的实数编码 |
| EliteCount | 2~4 | 保留每代最优个体，防止最优解丢失 |
| TolFun | 1e-6 | 适应度函数变化低于此值停止迭代 |
CrossoverFraction不是越大越好。太高意味着变异个体太少，算法容易集中到局部区域；太低则搜索效率下降。对于 5 变量问题，0.8 左右比较合适。MutationFcn我一般选用@mutationadaptfeasible，它会在边界约束内自适应调整变异步长，效果好于默认的高斯变异。EliteCount如果设为 0，最优解可能在迭代中被交叉破坏，导致收敛曲线回跳。
3.3 从单目标到多目标：Pareto前沿与折中解
单目标遗传算法只能给出一组最优解，但工程上往往需要在"响应快"和"控制能量小"之间找折中。比如汽车主动悬架，既要车身加速度小，又要悬架动挠度在限位范围内，两个指标互相冲突，这时单目标把两个指标加权成一个函数，权系数本身又成了需要调的东西。更直接的做法是用gamultiobj同时优化两个目标：第一个目标取超调量加调节时间，第二个目标取控制能量积分。
fitness_bi = @(vars) [lqr_time_score(vars, sys), lqr_energy_score(vars, sys)]; options_mo = optimoptions('gamultiobj', ... 'PopulationSize', 80, ... 'MaxGenerations', 100, ... 'ParetoFraction', 0.35); [vars_pareto, fvals] = gamultiobj(fitness_bi, nvars, [], [], [], [], lb, ub, [], options_mo);ParetoFraction控制返回的非支配解比例，0.35 表示最终种群里有 35% 的个体保留为 Pareto 前沿，数量一般在 20~30 个。从fvals里画出第一个目标对第二个目标的散点图，就能看到一条比较清晰的 Pareto 前沿，然后根据实际执行机构的功率上限选一个拐角点。注意gamultiobj默认使用拥挤距离排序，对多目标问题不需要额外归一化，但变量边界仍然要设，不然它会去探索毫无意义的负权重区域。
4. 优化过程中的三个坑与对策
4.1 可控性与LQR求解失败的保护
遗传算法随机产生的Q、R组合有很大概率不满足 LQR 求解前提。lqr要求A - B*K稳定，并且Q半正定、R正定。如果种群个体碰巧把某个q_i搜索到接近 0，可能导致黎卡提方程无解，lqr直接抛错。适应度函数里已经用try-catch挡了一层，但更早地判断可控性可以减少无效计算。
if rank(ctrb(sys.A, sys.B)) < size(sys.A, 1) score = 1e8; return; endctrb计算可控性矩阵，秩等于状态维数时才说明系统完全可控。注意这里检查的不是开环系统是否稳定，而是是否可控，因为 LQR 本身可以把开环不稳定系统镇定。对于可控性矩阵接近奇异的病态系统，用rank判断时还要设定一个容差，比如rank(ctrb(...), 1e-8)，否则数值误差可能导致误判。
4.1.1 不要用默认的tol判秩
Matlab 的rank默认容差是max(size(A))*eps(norm(A))量级，对系统矩阵本身数值差异大的实际问题太严格。比如状态量有角度和位移，量级差 100 倍，可控性矩阵的条件数就会变大，默认容差可能报告"不可控"，而实际上模型是可镇的。建议写成rank(ctrb(A,B), 1e-6)，把容差放宽到 1e-6，过滤掉真正的数值秩亏。
4.2 适应度函数里的控制量饱和
LQR 是基于线性模型的，它计算出的u = -Kx不感知执行机构的物理限幅。适应度函数如果只按线性仿真评估，遗传算法会倾向于给出很大的Q元素，让增益K变得很大，仿真曲线看起来收敛飞快，但实际u远超出执行机构能力。这属于"仿真看着最优，实物直接抖动"的典型问题。
u_sat = 5; u = -K * x; u = max(-u_sat, min(u_sat, u));在递推仿真里，每个时间步计算完u之后立刻加一个限幅钳位。对u做饱和处理后，系统的等效控制量不再是线性的，闭环性能会下降，遗传算法为了降低 ITAE，就得主动收敛到一组更温和的Q、R。这样就避免了理论增益过大问题。
4.2.1 饱和惩罚要加在什么位置
前面适应度函数里的u_norm是直接累加u^2，如果饱和钳位已经生效，累加时应该用饱和后的u。有些实现是先算线性u，再单独加一个"超出限幅的量"作为惩罚，其实不如统一使用饱和后的u自然，因为遗传算法优化的是实际物理过程，不是虚构的线性过程。
4.3 随机种子与重复性
遗传算法本身是随机算法，不同次运行得到的结果可能有明显差异。控制优化场景里，我们要的是可复现的设计，所以需要在主程序开头固定随机种子。但固定种子也有副作用：每次都一样，无法评估算法稳定性。
rng(42); % 固定种子得到确定结果 % rng('shuffle'); % 随机种子用于批量实验批量实验时，用for i = 1:10循环，每次用rng(i)，跑完 10 次后看最优 ITAE 的分布。如果 10 次最优值的最大值和最小值相差超过 20%，说明种群大小或者迭代代数不够，或者变量边界太宽，需要增大规模，而不是直接取某一次的结果。这是判断遗传算法收敛质量最廉价的实验方法。
5. 一个能直接用起来的优化模板
结合前面的段落，整理出下面这套模板，结构上只有一个主程序和两个适应度函数文件。拿到任何A、B矩阵，先改主程序前三行，再跑。
% 主程序模板 main_lqr_ga.m rng(42); A = ...; B = ...; sys = ss(A, B, eye(size(A,1)), 0); x0 = zeros(size(A,1), 1); x0(1) = 0.1; nvars = size(A,1) + size(B,2); lb = [-4*ones(1,size(A,1)), -6*ones(1,size(B,2))]; ub = [4*ones(1,size(A,1)), 2*ones(1,size(B,2))]; options = optimoptions('ga', 'PopulationSize', 60, 'MaxGenerations', 80); fitness = @(vars) lqr_fitness(vars, sys, x0, 1e-3, 5); [vars_opt, fval] = ga(fitness, nvars, [], [], [], [], lb, ub, [], options);模板的关键参数是x0的选取。x0代表系统的初始扰动，遗传算法是在这个特定扰动下寻找最优Q、R，换一个更大的初始扰动，最优解会变化。所以模板里的dt=1e-3、tfinal=5、x0(1)=0.1都不是随便设的，它们对应一个典型工况。实际使用时，应该把x0设成你最关心的工况，比如倒立摆的初始倾角 0.2 rad，或者悬架受到的路面冲击幅值。
另外一个实用技巧是把遗传算法的结果导入到 Simulink 里做非线性模型验证。Simulink 里把K_opt设为常量，用饱和模块和实际执行机构模型，仿真时间和模板里的tfinal一致，如果非线性模型下的 ITAE 和模板仿真差别超过 30%，说明线性化偏差太大，需要在适应度函数里加入工作点附近多个扰动点同时求平均，避免算法只盯着一个线性化点。这个验证步骤能直接暴露模板误差，也是这类优化设计工程里最容易被跳过的一环。先把模板跑通，再逐步加入饱和、扰动和第二个优化目标。
本文还有配套的精品资源，点击获取
