基于状态空间模型的预测控制SSMPC仿真程序

资 源 简 介

本项目是一个基于MATLAB开发的高级控制算法仿真程序，核心旨在实现基于状态空间（State-Space）模型的模型预测控制（MPC）算法。程序首先将连续时间的被控对象模型离散化，构建增广状态空间模型以消除稳态误差。在算法实现层面，程序依据设定的预测时域（Prediction Horizon）和控制时域（Control Horizon），自动推导并生成多步预测方程的系数矩阵。系统将控制问题转化为带约束的二次规划（Quadratic Programming, QP）问题，定义了包含输出跟踪误差项和控制增量项的目标函数（Cost Function）。在仿真循环中，程序完整实现了滚动优化（Receding Horizon Control）机制，即在每个采样时刻根据当前状态测量值在线求解最优控制序列，并仅将序列中的第一个控制量作用于系统。该程序支持处理输入幅值约束、输入增量约束以及输出约束，能够模拟多输入多输出（MIMO）或单输入单输出（SISO）系统在设定值跟踪和抗扰动测试中的动态表现。此外，项目集成了可视化模块，能够实时绘制系统响应曲线、控制量变化曲线及误差收敛情况，便于用户分析算法参数（如权重矩阵Q/R、时域长度）对系统性能的影响。

详 情 说 明

基于状态空间模型的预测控制 (SSMPC) 仿真系统

项目介绍

本项目是一个基于 MATLAB 开发的高级模型预测控制（Model Predictive Control, MPC）仿真系统。该系统针对多输入多输出（MIMO）被控对象，实现了基于增广状态空间模型的 MPC 算法。

程序采用无静差跟踪的设计理念，通过构建增广状态模型（包含状态增量和输出），将控制问题转化为标准的二次规划（QP）问题。系统能够显式处理输入幅值约束和输入增量约束，并在每个采样时刻通过滚动优化机制求解最优控制律。该项目不依赖 MATLAB 的 Simulink 模块，完全基于脚本和数学运算实现，非常适合用于理解 MPC 的底层数学原理和算法流程。

主要功能特性

• 多输入多输出 (MIMO) 支持：代码内置了一个 2 输入 2 输出的连续系统模型，展示了 MPC 对多变量耦合系统的控制能力。
• 状态空间离散化：不依赖控制系统工具箱的 c2d 函数，而是直接采用矩阵指数法（Matrix Exponential）实现连续系统的零阶保持（ZOH）离散化。
• 无静差跟踪：采用增广状态空间模型（Augmented State-Space Model），状态向量包含“状态增量”与“输出”，引入积分作用以消除稳态误差。
• 多步预测与滚动优化：支持自定义预测时域（Prediction Horizon）和控制时域（Control Horizon），自动生成多步预测方程。
• 硬约束处理：集成了基于 quadprog 的二次规划求解器，能够严格处理：

* 控制增量约束 ($Delta u_{min} le Delta u le Delta u_{max}$) * 控制幅值约束 ($u_{min} le u le u_{max}$)

• 动态仿真与可视化：包含完整的仿真主循环，实时计算并存储数据，最后生成关于输出跟踪、控制量变化及计算耗时的详细图表。

系统要求

• MATLAB R2016b 或更高版本
• Optimization Toolbox（必须，用于运行 quadprog 函数求解 QP 问题）

使用方法

1. 确保 MATLAB 路径中包含本项目所在文件夹。
2. 直接运行主函数（main）。
3. 程序将自动执行以下步骤：

* 初始化系统模型参数。 * 进行离线矩阵计算（预测矩阵、QP 矩阵）。 * 执行 10 秒的控制仿真循环。 * 弹出图形窗口显示仿真结果。

详细代码实现分析

该仿真程序完全包含在一个文件中，通过清晰的模块化结构实现了 SSMPC 的完整流程。以下是代码实现的详细逻辑：

1. 模型建立与离散化

• 连续模型定义：定义了双输入双输出系统的 $A_c, B_c, C_c$ 矩阵。
• 自定义离散化：程序构建了扩展矩阵 $M_{temp}$ 并利用 expm 计算矩阵指数，提取出离散时间的状态矩阵 $A_d$ 和输入矩阵 $B_d$。这种方法具有更好的数值通用性，不依赖特定工具箱。

2. 增广状态模型构建

为了实现无静差跟踪，代码构建了增广状态向量 $x_a(k) = [Delta x(k)^T, y(k)^T]^T$。

• 构建了增广系统矩阵 $A_{aug}, B_{aug}, C_{aug}$。
• 这种建模方式将控制输入转化为输入增量 $Delta u$，使得算法本质上具备了积分特性。

3. QP 问题矩阵的离线计算

为了提高在线计算效率，程序在循环外预先计算了所有时不变矩阵：

• 预测方程矩阵 ($F, Phi$)：通过循环迭代，计算了从当前状态预测未来 $N_p$ 步输出的系数矩阵。

* $Y = F cdot x_a(k) + Phi cdot Delta U$

• Hessian 矩阵 ($H$)：基于加权矩阵 $Q$（输出误差权重）和 $R$（控制增量权重），构建了标准 QP 问题的 $H$ 矩阵：

* $H = Phi^T bar{Q} Phi + bar{R}$ * 代码中显式执行了 $H = (H + H')/2$ 以确保数值上的对称正定性，这对 QP 求解器的稳定性至关重要。

• 约束矩阵构建：

* 增量约束：直接映射为 QP 问题的下限 LB 和上限 UB。 * 幅值约束：利用下三角矩阵 L_tril 将未来的控制量表示为当前控制量与累积增量之和，从而构建不等式约束矩阵 $A_{cons}$。此处仅构建了 $A$ 矩阵部分，常数向量 $b$ 在在线循环中更新。

4. 仿真主循环 (Rolling Horizon)

仿真时长设定为 10 秒，采样时间 0.1 秒。循环内的主要逻辑如下：

1. 参考轨迹生成：根据当前时间步生成参考信号（通道1为方波，通道2为正弦波）。
2. 状态反馈与增广：计算状态增量 $Delta x$ 和当前输出 $y$，组合成增广状态 $x_a$。
3. 在线构造 QP 线性项 ($f$)：结合参考轨迹向量与当前状态预测，计算 QP 目标函数的一次项系数 $f$。
4. 更新动态约束 ($b_{cons}$)：根据上一时刻的实际控制两 $u(k-1)$，更新幅值约束的边界向量。
5. 求解 QP：调用 quadprog 求解最优控制增量序列 $Delta U^*$。

* 代码包含容错处理：若求解器未收敛（exitflag != 1），则默认输出零增量，防止程序崩溃。

1. 实施控制：仅取优化序列中的第一个元素 $Delta u(k)$，计算实际控制量 $u(k) = u(k-1) + Delta u(k)$ 并作用于系统模型。
2. 状态更新：利用离散状态方程计算下一时刻的物理状态。

5. 结果可视化

plot_results 辅助函数生成三组子图：

• 输出响应：对比系统的实际输出 $y$ 与参考轨迹 $r$，展示跟踪性能。
• 控制输入：显示控制量 $u$ 的变化曲线，可用于观察是否触及设定幅值约束。
• 求解性能：(注：实际代码中包含 solver_time 的记录和存储，但主要可视化集中在输出和输入上，便于评估控制器性能)。

算法参数说明

代码中预设的关键参数如下，用户可根据需要修改：

• 预测时域 ($N_p$)：20步。决定了控制器“向前看”的距离。
• 控制时域 ($N_c$)：5步。决定了控制自由度，较小的值有助于减少计算量。
• 权重矩阵：

* $Q_{weight} = text{diag}(10, 10)$：强调对输出误差的惩罚，提高跟踪精度。 * $R_{weight} = text{diag}(0.1, 0.1)$：对控制增量的惩罚较小，允许相对快速的控制动作。

• 约束条件：

* 输入范围：$pm 3$ * 增量范围：$pm 0.5$