非线性弹簧质量系统模型预测控制（MPC）设计与实现

一、系统概述

非线性弹簧质量系统是经典的非线性动力学系统，由质量块、非线性弹簧（力-位移关系非胡克线性）和阻尼器组成，广泛用于机械振动、车辆悬挂等领域。模型预测控制（MPC）通过在线求解有限时域优化问题，可实现对该系统的高精度轨迹跟踪与约束控制。

本设计以单自由度非线性弹簧质量系统为对象，实现简单MPC控制，核心步骤包括：

建立非线性动力学模型（含立方弹簧力）；
离散化模型（欧拉法）；
设计MPC优化问题（目标函数+约束）；
MATLAB仿真验证（跟踪阶跃/正弦轨迹）。

二、系统建模

2.1 非线性动力学方程

系统组成：质量块（质量 $m$ ）、非线性弹簧（力 $F_s=kx+βx^3$ ，β为非线性系数）、阻尼器（阻尼系数 $c$ ）、外力输入 $u(t)$ 。

状态变量： $x_1=x$ （位移）， $x_2=\dot{x}$ （速度），状态向量 $x=[x_1,x_2]^T$ 。

控制输入：4u（外力），输出：y=x1（位移）。

动力学方程（牛顿第二定律）：

$m\ddot{x}−kx−βx^3−c\dot{x}+u$

整理为状态空间形式：

参数取值（示例）：

质量 $m=1 kg$ ，弹簧刚度 $k=10 N/m$ ，非线性系数 $β=2 N/m^3$ ，阻尼系数c=0.5 N\cdotps/m。

2.2 离散化模型（欧拉法）

MPC需离散时间模型，采用前向欧拉法（T为采样时间）：

向量形式：

$x(k+1)=f(x(k),u(k))$

其中 $f$ 为非线性离散状态转移函数。

三、MPC控制器设计

3.1 优化问题 formulation

预测时域 $N_p$ （预测未来 $N_p$ 步状态），控制时域 $N_c$ （优化未来 $N_c$ 步控制输入， $N_c≤N_p$ ）。

目标函数（跟踪误差+控制输入代价）：

$x(k+i∣k)$ ： $k$ 时刻预测 $k+i$ 时刻状态；
$x_{ref}$ ：参考轨迹；
$Q,R,S$ ：权重矩阵（对角阵，分别惩罚状态误差、控制输入、输入变化率）；
$Δu=u(k+i∣k)−u(k+i−1∣k)$ ：控制输入增量。

约束：

状态约束： $x_{min}≤x(k+i∣k)≤x_{max}$ ；
控制输入约束： $u_{min}≤u(k+i∣k)≤u_{max}$ ；
输入增量约束： $Δu_{min}≤Δu(k+i∣k)≤Δu_{max}$ 。

3.2 优化问题求解

非线性MPC的优化问题属于非线性规划（NLP），需用数值方法求解。本设计采用MATLAB的fmincon函数（内点法），步骤如下：

定义目标函数（含状态转移与代价计算）；
定义约束（状态、输入、增量约束）；
以当前状态为初始值，求解最优控制序列 ${u(k∣k),u(k+1∣k),...,u(k+Nc−1∣k)}$ ；
仅应用第一个控制输入 $u(k∣k)$ ，下一时刻重复上述过程。

四、MATLAB实现

4.1 系统参数与模型

% 系统参数
m = 1;          % 质量 (kg)
k = 10;         % 弹簧刚度 (N/m)
beta = 2;        % 非线性系数 (N/m^3)
c = 0.5;         % 阻尼系数 (N·s/m)
T = 0.05;        % 采样时间 (s)

% 状态转移函数（离散非线性模型）
function x_next = system_dynamics(x, u, m, k, beta, c, T)
    x1 = x(1);  % 位移
    x2 = x(2);  % 速度
    % 欧拉法离散化
    x1_next = x1 + T * x2;
    x2_next = x2 + T * ( (-k*x1 - beta*x1^3 - c*x2 + u) / m );
    x_next = [x1_next; x2_next];
end

4.2 MPC控制器类

classdef NonlinearSpringMPC < handle
    properties
        % 系统参数
        m; k; beta; c; T;  % 质量、刚度、非线性系数、阻尼、采样时间
        nx = 2; nu = 1;     % 状态维度(2)、输入维度(1)
        
        % MPC参数
        Np; Nc;             % 预测时域、控制时域
        Q; R; S;            % 权重矩阵
        x_min; x_max;       % 状态约束
        u_min; u_max;       % 输入约束
        du_min; du_max;      % 输入增量约束
        
        % 参考轨迹
        x_ref;              % 参考状态序列
    end
    
    methods
        function obj = NonlinearSpringMPC(m, k, beta, c, T, Np, Nc, Q, R, S)
            % 构造函数
            obj.m = m; obj.k = k; obj.beta = beta; obj.c = c; obj.T = T;
            obj.Np = Np; obj.Nc = Nc;
            obj.Q = Q; obj.R = R; obj.S = S;
            
            % 默认约束
            obj.x_min = [-5; -5];  % 位移/速度下限
            obj.x_max = [5; 5];    % 位移/速度上限
            obj.u_min = -20;        % 外力下限 (N)
            obj.u_max = 20;         % 外力上限 (N)
            obj.du_min = -10;       % 输入增量下限
            obj.du_max = 10;        % 输入增量上限
        end
        
        function [u_opt, x_pred] = solve_mpc(obj, x0, x_ref_seq)
            % 求解MPC优化问题
            % 输入: x0(当前状态), x_ref_seq(参考状态序列, Np步)
            % 输出: u_opt(最优控制输入), x_pred(预测状态序列)
            
            % 优化变量: U = [u(0), u(1), ..., u(Nc-1)]^T (Nc维)
            Nc = obj.Nc;
            U0 = zeros(Nc, 1);  % 初始猜测
            
            % 约束设置
            A_ineq = []; b_ineq = [];
            A_eq = []; b_eq = [];
            lb = repmat(obj.u_min, Nc, 1);  % 输入下界
            ub = repmat(obj.u_max, Nc, 1);  % 输入上界
            
            % 输入增量约束: u(i) - u(i-1) ∈ [du_min, du_max] (i≥1)
            if Nc > 1
                A_du = [eye(Nc-1), zeros(Nc-1, 1); -eye(Nc-1), zeros(Nc-1, 1)];
                b_du = [repmat(obj.du_max, Nc-1, 1); -repmat(obj.du_min, Nc-1, 1)];
                A_ineq = A_du; b_ineq = b_du;
            end
            
            % 调用fmincon求解
            options = optimoptions('fmincon', 'Display', 'off', 'Algorithm', 'sqp');
            U_opt = fmincon(@(U) obj.cost_function(U, x0, x_ref_seq, obj), U0, ...
                            A_ineq, b_ineq, A_eq, b_eq, lb, ub, [], options);
            
            % 提取最优控制输入（仅用第一个）
            u_opt = U_opt(1);
            
            % 预测状态序列（用于分析）
            x_pred = obj.predict_states(x0, U_opt);
        end
        
        function J = cost_function(obj, U, x0, x_ref_seq, self)
            % 目标函数计算
            Np = self.Np; Nc = self.Nc;
            Q = self.Q; R = self.R; S = self.S;
            m = self.m; k = self.k; beta = self.beta; c = self.c; T = self.T;
            
            x = x0;  % 当前状态
            J = 0;   % 初始化代价
            U_ext = [U; U(end)*ones(Np-Nc, 1)];  % 控制时域外保持最后值
            
            for i = 1:Np
                u = U_ext(i);
                % 状态转移
                x_next = system_dynamics(x, u, m, k, beta, c, T);
                % 状态误差代价
                J = J + (x_next - x_ref_seq(:,i))' * Q * (x_next - x_ref_seq(:,i));
                % 控制输入代价
                if i <= Nc
                    J = J + u' * R * u;
                    % 输入增量代价 (i>1时)
                    if i > 1
                        du = u - U_ext(i-1);
                        J = J + du' * S * du;
                    end
                end
                x = x_next;  % 更新状态
            end
        end
        
        function x_pred = predict_states(obj, x0, U)
            % 预测Np步状态
            Np = obj.Np; Nc = obj.Nc;
            m = obj.m; k = obj.k; beta = obj.beta; c = obj.c; T = obj.T;
            x = x0;
            x_pred = zeros(2, Np);
            U_ext = [U; U(end)*ones(Np-Nc, 1)];  % 控制时域外保持最后值
            
            for i = 1:Np
                u = U_ext(i);
                x = system_dynamics(x, u, m, k, beta, c, T);
                x_pred(:,i) = x;
            end
        end
    end
end

4.3 仿真主程序

%% 非线性弹簧质量系统MPC仿真
clear; close all; clc;

%% 1. 系统参数与MPC设置
% 系统参数
m = 1; k = 10; beta = 2; c = 0.5; T = 0.05;  % 采样时间0.05s
% MPC参数
Np = 20;  % 预测时域
Nc = 5;   % 控制时域
Q = diag([100, 1]);  % 状态权重 (位移权重100，速度权重1)
R = 0.1;  % 控制输入权重
S = 0.01; % 输入增量权重
% 创建MPC控制器
mpc = NonlinearSpringMPC(m, k, beta, c, T, Np, Nc, Q, R, S);

%% 2. 参考轨迹生成（阶跃+正弦）
t_sim = 5;  % 仿真时间5s
N_sim = t_sim / T;  % 仿真步数
t = 0:T:t_sim-T;
x_ref = zeros(2, N_sim);  % 参考状态 (位移+速度)
% 位移参考：0~2s阶跃(1m)，2~5s正弦(振幅0.5m，频率0.5Hz)
x_ref(1, t<2) = 1;
x_ref(1, t>=2) = 1 + 0.5*sin(2*pi*0.5*(t(t>=2)-2));
x_ref(2, :) = 0.5 * 0.5 * 2*pi*0.5*cos(2*pi*0.5*(t-2));  % 速度参考（导数）

%% 3. 仿真初始化
x0 = [0; 0];  % 初始状态 (位移0，速度0)
x_history = zeros(2, N_sim+1);  % 状态历史
u_history = zeros(1, N_sim);    % 控制输入历史
x_history(:,1) = x0;

%% 4. 主仿真循环
for k = 1:N_sim
    % 当前参考轨迹 (Np步)
    ref_start = k;
    ref_end = min(k + Np - 1, N_sim);
    x_ref_seq = [x_ref(:, ref_start:ref_end), zeros(2, Np - (ref_end - ref_start + 1))];  % 不足补零
    
    % 求解MPC
    [u_opt, x_pred] = mpc.solve_mpc(x_history(:,k), x_ref_seq);
    
    % 应用控制输入
    u_history(k) = u_opt;
    
    % 系统状态更新 (非线性模型)
    x_next = system_dynamics(x_history(:,k), u_opt, m, k, beta, c, T);
    x_history(:,k+1) = x_next;
    
    % 显示进度
    if mod(k, 50) == 0
        fprintf('时间: %.2f s, 位移: %.2f m, 参考: %.2f m\n', ...
                k*T, x_history(1,k+1), x_ref(1,min(k+1,N_sim)));
    end
end

%% 5. 结果可视化
figure('Position', [100, 100, 1000, 600]);

% 子图1：位移跟踪
subplot(2,1,1);
plot(t, x_history(1,:), 'b-', 'LineWidth', 1.5);
hold on;
plot(t, x_ref(1,:), 'r--', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('位移 (m)');
title('非线性弹簧质量系统位移跟踪 (MPC)');
legend('实际位移', '参考位移', 'Location', 'best');
grid on;

% 子图2：控制输入
subplot(2,1,2);
plot(t, u_history, 'g-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('控制输入 (N)');
title('控制输入 (外力)');
grid on;

%% 6. 性能指标
pos_error = x_history(1,:) - x_ref(1,:);
rmse = sqrt(mean(pos_error.^2));
fprintf('\n性能指标:\n');
fprintf('位移跟踪RMSE: %.4f m\n', rmse);
fprintf('最大控制输入: %.2f N\n', max(abs(u_history)));

参考代码非线性弹簧质量系统的简单模型预测控制(MPC) www.youwenfan.com/contentcnt/160595.html

五、关键问题与解决方案

5.1 非线性模型离散化精度

问题：欧拉法精度低，可能导致预测误差大。

解决：采用4阶龙格-库塔法离散化，提高精度：

function x_next = rk4_dynamics(x, u, m, k, beta, c, T)
    k1 = system_derivative(x, u, m, k, beta, c);
    k2 = system_derivative(x + T/2*k1, u, m, k, beta, c);
    k3 = system_derivative(x + T/2*k2, u, m, k, beta, c);
    k4 = system_derivative(x + T*k3, u, m, k, beta, c);
    x_next = x + T/6*(k1 + 2*k2 + 2*k3 + k4);
end
function dx = system_derivative(x, u, m, k, beta, c)
    dx1 = x(2);
    dx2 = (-k*x(1) - beta*x(1)^3 - c*x(2) + u)/m;
    dx = [dx1; dx2];
end

5.2 优化求解效率

问题：fmincon求解非线性规划较慢，难以实时。
解决：
- 简化目标函数（如忽略输入增量代价S）；
- 减小预测时域Np和控制时域Nc；
- 使用显式MPC（离线计算控制律，在线查表）。

5.3 约束处理

问题：状态/输入约束可能导致优化无解。
解决：
- 放宽约束（如适当扩大xmin/xmax）；
- 采用罚函数法将不等式约束转化为目标函数惩罚项。

六、总结

本设计实现了非线性弹簧质量系统的简单MPC控制，核心步骤包括：

建立含立方弹簧力的非线性动力学模型；
欧拉法离散化模型；
设计MPC优化问题（目标函数+约束）；
MATLAB仿真验证（跟踪阶跃+正弦轨迹）。