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

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

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

LQR 控制器基本原理

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

$$J = \int (x^T Q x + u^T R u) dt$$

其中：

• $Q \ge 0$：状态加权矩阵
• $R > 0$：控制输入加权矩阵
• $u = -Kx$：状态反馈控制律

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

$$K = R^{-1} B^T P$$ $$A P + P A^T - P B R^{-1} B^T P + Q = 0$$

遗传算法优化 LQR 设计

优化问题分析

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

算法流程

1. 初始化种群
2. 评估适应度
3. 判断是否满足终止条件？
   • 是：输出最优解
   • 否：执行选择操作 -> 交叉操作 -> 变异操作 -> 生成新一代种群 -> 返回步骤 2

MATLAB 实现

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

function [fitness, stability] = lqr_fitness_function(params, A, B, Q_structure, R_structure)
    % 提取参数
    n = size(A, 1);
    m = size(B, 2);
    
    % 重构 Q 矩阵
    if strcmp(Q_structure, 'diagonal')
        q_diag = params(1:n);
        Q = diag(q_diag);
        param_offset = n;
    elseif strcmp(Q_structure, 'full')
        Q = reshape(params(1:n*n), n, n);
        Q = (Q + Q') / 2; % 确保对称性
        param_offset = n*n;
    end
    
    % 重构 R 矩阵
    if strcmp(R_structure, 'scalar')
        R = params(param_offset + 1);
        R = R * eye(m);
    elseif strcmp(R_structure, 'diagonal')
        r_diag = params(param_offset + 1:param_offset + m);
        R = diag(r_diag);
        param_offset = param_offset + m;
    elseif strcmp(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;
    end
end

function 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');
    if isempty(idx)
        ts = time(end);
    else
        ts = time(idx);
    end
end

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); % 输入维度
    
    % 根据矩阵结构确定变量个数
    if strcmp(Q_structure, 'diagonal')
        nQ = n;
    elseif strcmp(Q_structure, 'full')
        nQ = n*n;
    end
    
    if strcmp(R_structure, 'scalar')
        nR = 1;
    elseif strcmp(R_structure, 'diagonal')
        nR = m;
    elseif strcmp(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);
        parfor i = 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));
        for i = 1:pop_size
            candidates = randperm(pop_size, 2);
            [~, winner] = min(fitness(candidates));
            new_population(i,:) = population(candidates(winner), :);
        end
        
        % 交叉操作 (算术交叉)
        for i = 1:2:pop_size-1
            if 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;
            end
        end
        
        % 变异操作 (高斯变异)
        for i = 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);
            end
        end
        population = new_population;
    end
    
    % 提取最优解
    best_params = best_individual;
    [~, stability] = lqr_fitness_function(best_params, A, B, Q_structure, R_structure);
    
    % 重构最优 Q 和 R
    if strcmp(Q_structure, 'diagonal')
        q_diag = best_params(1:n);
        Q_opt = diag(q_diag);
        param_offset = n;
    elseif strcmp(Q_structure, 'full')
        Q_opt = reshape(best_params(1:n*n), n, n);
        Q_opt = (Q_opt + Q_opt') / 2;
        param_offset = n*n;
    end
    
    if strcmp(R_structure, 'scalar')
        R_opt = best_params(param_offset + 1) * eye(m);
    elseif strcmp(R_structure, 'diagonal')
        r_diag = best_params(param_offset + 1:param_offset + m);
        R_opt = diag(r_diag);
    elseif strcmp(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 = [0 1 0 0; 0 0 -(m*g)/M 0; 0 0 0 1; 0 0 ((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);
if strcmp(Q_struct, 'diagonal')
    q_diag = best_params(1:n);
    Q_opt = diag(q_diag);
    R_opt = best_params(n+1);
elseif strcmp(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;
    end
end

性能评估与可视化

function plot_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;
    for i = 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));
    
    for i = 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

应用场景与优势

典型应用领域

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

算法优势

特性	传统 LQR	GA-LQR
全局最优性	局部最优	全局搜索
约束处理	困难	灵活
非线性系统	不适用	可扩展
参数调整	经验依赖	自动优化
计算复杂度	低	较高

实际应用建议

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