简介这是一份基于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(n1)/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权重对数-44对应权重 1e-4 到 1e4lg(q2)状态2权重对数-44同上lg(r1)控制权重对数-62控制代价通常比状态权重小边界不一定要对称如果控制量特别珍贵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(n1: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(n1: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)固定随机种子保证每次运行结果一致否则遗传算法的初始种群是随机的同样参数可能跑出不同结果不利于对比实验。nvars5对应四个状态权重和一个控制权重。UseParallel设为true可以在多核 CPU 上并行计算适应度但需要提前用parpool开启并行池否则这个选项会被忽略。这里被控对象的A、B是倒立摆在平衡点的线性化结果换成自己的模型时只需要改A、B和初值。3.2 关键参数配置种群、交叉变异、边界约束遗传算法的效果很大程度卡在参数选择上。ga的默认参数对控制器优化这类中等规模问题并不理想建议按下面的表调整参数名建议值作用PopulationSize40~80太小容易早熟太大收敛慢MaxGenerations60~120每一代耗时和收敛判断的平衡CrossoverFraction0.7~0.9交叉比例高种群多样性好MutationFcnmutationadaptfeasible自适应变异适合带边界约束的实数编码EliteCount2~4保留每代最优个体防止最优解丢失TolFun1e-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换一个更大的初始扰动最优解会变化。所以模板里的dt1e-3、tfinal5、x0(1)0.1都不是随便设的它们对应一个典型工况。实际使用时应该把x0设成你最关心的工况比如倒立摆的初始倾角 0.2 rad或者悬架受到的路面冲击幅值。另外一个实用技巧是把遗传算法的结果导入到 Simulink 里做非线性模型验证。Simulink 里把K_opt设为常量用饱和模块和实际执行机构模型仿真时间和模板里的tfinal一致如果非线性模型下的 ITAE 和模板仿真差别超过 30%说明线性化偏差太大需要在适应度函数里加入工作点附近多个扰动点同时求平均避免算法只盯着一个线性化点。这个验证步骤能直接暴露模板误差也是这类优化设计工程里最容易被跳过的一环。先把模板跑通再逐步加入饱和、扰动和第二个优化目标。本文还有配套的精品资源点击获取