如何在Matlab仿真中为波音747预测控制模型分配约束值及确定上下限？

问题背景 #

本案例针对40000英尺高度飞行的波音747开发预测控制模型，通过调整发动机推力与升降舵，维持最优空速和爬升率。由于输出（空速、爬升率）采用相对于标称参考值的偏差变量（而非绝对值）评估，输入（推力、升降舵）的约束也需遵循相同的偏差逻辑。核心需求是：在Matlab仿真中为输入分配与输出一致的偏差约束，并确定合理的上下限。

模型基础 #

该模型基于波音747在40000英尺、0.8马赫（774ft/s）水平飞行状态下的线性化状态空间模型：

• 状态变量：水平速度偏差、垂直速度偏差、角速度、俯仰角（单位：crad=0.01rad）
• 输入变量：升降舵偏差、推力偏差
• 输出变量：空速偏差、爬升率偏差

Gekko示例代码（Python） #

以下是实现该预测控制的Gekko参考代码：

from gekko import GEKKO
import numpy as np

## 波音747线性模型
# 飞行状态：40000英尺高度水平飞行，速度774ft/s（0.8马赫）
# 状态量：u-uw(ft/s)水平速度偏差、w-ww(ft/s)垂直速度偏差、q(crad/s)角速度、theta(crad)俯仰角
# 输入量：e升降舵偏差、t推力偏差
# 输出量：空速偏差(u-uw)、爬升率偏差(-w + 774*theta)

A = np.array([[-.003, 0.039, 0, -0.322],
              [-0.065, -0.319, 7.74, 0],
              [0.020, -0.101, -0.429, 0],
              [0, 0, 1, 0]])

B = np.array([[0.01, 1],
              [-0.18, -0.04],
              [-1.16, 0.598],
              [0, 0]])

C = np.array([[1, 0, 0, 0],
              [0, -1, 0, 7.74]])

# 构建Gekko状态空间模型
m = GEKKO()
x,y,u = m.state_space(A,B,C,D=None)

m.time = [0, 0.1, 0.2, 0.4, 1, 1.5, 2, 3, 4, 5, 6, 7, 8, 10, 12, 15]
m.options.imode = 6  # 预测控制模式
m.options.nodes = 3

## 操纵变量(MV)约束与调优
# 为升降舵、推力偏差设置上下限、移动成本
for i in range(len(u)):
    u[i].lower = -5    # 偏差下限
    u[i].upper = 5     # 偏差上限
    u[i].dcost = 1     # 移动惩罚系数
    u[i].status = 1    # 启用该操纵变量

## 受控变量(CV)调优
# 轨迹跟踪时间常数
y[0].tau = 5    # 空速偏差轨迹时间常数
y[1].tau = 8    # 爬升率偏差轨迹时间常数
# 轨迹初始化模式：2=随周期重新对齐中心
y[0].tr_init = 2
y[1].tr_init = 2
# 设定偏差目标范围（死区）
y[0].sphi= -8.5
y[0].splo= -9.5
y[1].sphi= 5.4
y[1].splo= 4.6

y[0].status = 1
y[1].status = 1

m.solve()

# 读取轨迹结果
import json
with open(m.path+'//results.json') as f:
    results = json.load(f)
air_speed = y[0].name
climb_rate = y[1].name

# 绘图
import matplotlib.pyplot as plt
plt.figure(1)
plt.subplot(311)
plt.plot(m.time,u[0],'b-',lw=2.0)
plt.plot(m.time,u[1],'g:',lw=2.0)
plt.legend(['升降舵偏差','推力偏差'])
plt.ylabel('操纵变量')

plt.subplot(312)
plt.plot(m.time,y[0],'r-',lw=2.0)
plt.plot(m.time,results[air_speed+'.tr_hi'],'k:')
plt.plot(m.time,results[air_speed+'.tr_lo'],'k:')
plt.legend(['空速偏差','轨迹上限','轨迹下限'])
plt.ylabel('空速偏差(ft/s)')

plt.subplot(313)
plt.plot(m.time,y[1],'r-',lw=2.0)
plt.plot(m.time,results[climb_rate+'.tr_hi'],'k:')
plt.plot(m.time,results[climb_rate+'.tr_lo'],'k:')
plt.legend(['爬升率偏差','轨迹上限','轨迹下限'])
plt.ylabel('爬升率偏差(ft/s)')

plt.show()

Matlab中输入偏差约束的实现方法 #

1. 核心思路：统一偏差变量逻辑 #

由于输出采用相对于标称参考值的偏差，输入必须同步采用偏差变量建模：

• 定义实际输入 = 标称参考输入 + 输入偏差变量
• 所有约束直接施加在输入偏差变量上，而非绝对输入值

2. 输入偏差上下限的设定步骤 #

步骤1：确定标称输入值

首先获取标称飞行状态（40000英尺、0.8马赫水平飞行）下的基准输入值：

• 升降舵标称值：水平飞行时的平衡偏转角度（通常为0或小角度，可从模型或飞机手册获取）
• 推力标称值：维持水平匀速飞行所需的推力（可通过模型稳态解计算，或参考飞机手册）

步骤2：推导偏差上下限

有两种方法确定偏差范围：

方法A：基于输入物理极限

1. 从飞机手册或模型特性获取输入的绝对物理范围：
   • 升降舵：例如物理偏转范围为[-15°, +15°]（需转换为模型单位，如crad）
   • 推力：例如推力百分比范围为[20%, 100%]（转换为模型对应的数值单位）
2. 计算偏差范围：

   输入偏差上限 = 输入绝对上限 - 标称输入值
   输入偏差下限 = 输入绝对下限 - 标称输入值
   

方法B：基于输出允许偏差反推

根据输出（空速、爬升率）的允许偏差范围，结合状态空间模型的B矩阵（输入对状态的影响）或输入输出灵敏度，反推输入偏差的合理范围，避免输出超出约束。例如：

• 若空速允许偏差为±10ft/s，通过模型计算升降舵/推力对空速的增益，反推输入偏差的最大允许值。

3. Matlab仿真中的具体实现 #

以Matlab MPC工具箱为例：

步骤1：构建偏差形式的状态空间模型

将输入、输出均定义为偏差变量，模型矩阵直接使用案例中的A/B/C（已为线性化偏差模型）：

A = [-.003, 0.039, 0, -0.322;
     -0.065, -0.319, 7.74, 0;
     0.020, -0.101, -0.429, 0;
     0, 0, 1, 0];
B = [0.01, 1;
     -0.18, -0.04;
     -1.16, 0.598;
     0, 0];
C = [1, 0, 0, 0;
     0, -1, 0, 7.74];
D = zeros(2,2);

mpcobj = mpc(ss(A,B,C,D));

步骤2：设置输入偏差约束

假设已计算出：

• 升降舵偏差上下限：[-5, 5]（对应物理范围相对于标称值的偏移）
• 推力偏差上下限：[-5, 5]
  直接为操纵变量设置约束：

% 设定升降舵（第1个MV）的偏差约束
mpcobj.MV(1).Min = -5;
mpcobj.MV(1).Max = 5;
% 设定推力（第2个MV）的偏差约束
mpcobj.MV(2).Min = -5;
mpcobj.MV(2).Max = 5;

% 可选：设置输入移动速率约束（避免动作过于剧烈）
mpcobj.MV(1).RateMin = -2;
mpcobj.MV(1).RateMax = 2;
mpcobj.MV(2).RateMin = -1;
mpcobj.MV(2).RateMax = 1;

步骤3：设置输出偏差目标

将输出参考值设为偏差目标（如空速偏差目标范围[-9.5, -8.5]，爬升率偏差目标范围[4.6, 5.4]）：

% 设置空速偏差（第1个CV）的死区目标
mpcobj.CV(1).Lower = -9.5;
mpcobj.CV(1).Upper = -8.5;
% 设置爬升率偏差（第2个CV）的死区目标
mpcobj.CV(2).Lower = 4.6;
mpcobj.CV(2).Upper = 5.4;

步骤4：运行仿真

T = 0:0.1:15;
r = [repmat([-9 5], length(T), 1)]; % 输出偏差目标值
[y,t,u] = sim(mpcobj, T, r);

% 绘图
figure;
subplot(3,1,1);
plot(t, u(:,1), 'b-', t, u(:,2), 'g:');
legend('升降舵偏差','推力偏差');
ylabel('操纵变量');

subplot(3,1,2);
plot(t, y(:,1), 'r-');
hold on;
plot(t, repmat(-9.5, length(t),1), 'k--', t, repmat(-8.5, length(t),1), 'k--');
legend('空速偏差','目标下限','目标上限');
ylabel('空速偏差(ft/s)');

subplot(3,1,3);
plot(t, y(:,2), 'r-');
hold on;
plot(t, repmat(4.6, length(t),1), 'k--', t, repmat(5.4, length(t),1), 'k--');
legend('爬升率偏差','目标下限','目标上限');
ylabel('爬升率偏差(ft/s)');

内容的提问来源于stack exchange，提问作者Andrea Emanuele Gambardella