基于预测控制的汽轮机转速调节

参考文献
汽轮机调节系统建模：李京,张志荣.汽轮机调节系统状态空间方程动态分析[J].热力发电,2002,31(6):22-25.
预测控制原理和代码：https://www.bilibili.com/video/BV1cL411n7KV
https://blog.csdn.net/qq_44940689/article/details/139808413

为了学习预测控制的仿真，以汽轮机调节系统为案例进行控制，所以在原代码基础上根据具体情况需要修改了。
对汽轮机调节系统并不了解，所以我的控制输入加入系统的地方纯属自找，自推理，自想想，方便控制，无依据。

目录

  1. 机理法建立状态空间模型
    1.1 汽轮机转速调节系统运动的微分方程
    1.1.1 汽轮发电机转子的运动方程
    1.1.2 蒸汽容积方程
    1.1.3 调速器方程
    1.1.4 油动机方程
    1.1.5 错油门运动方程
    1.2 汽轮机转速调节系统的状态空间模型

  2. 基于模型预测控制的汽轮机转速控制
    2.1. 最优控制问题
    2.2. 模型预测控制原理
    2.3. 基于汽轮机转速控制的MPC控制器实现
    2.3.1. 目标函数
    2.3.2. 二次规划型目标函数

  3. 算例仿真
    3.1. 参数设置与仿真程序
    3.1.1. MPC程序
    3.1.2. PID仿真模块
    3.2. MPC控制
    3.2.1. 输出权重和输入权重
    3.2.2. 其他状态权重
    3.3. PID控制与对比
    3.4. 模型不确定性下的MPC性能

  4. 结论

附录
附录1：MPC.m程序
附录2：MPC_Matrices.m程序
附录3：Prediction.m程序

基于预测控制的汽轮机转速调节器

1 机理法建立状态空间模型
利用机理法分析汽轮机调节系统的动态特性，首先必须对汽轮机的运行原理有相当的了解, 运用一些相应的物理基本定律对汽轮机系统的内部机理进行理论上的分析来推导出数学模型，也就是研究并写出系统的微分方程。而组成汽轮机调节系统的各个环节不论是机械的、液压的、电气的或是热力的都可以用微分方程加以描述。此时得到了多个微分方程，需要用状态空间模型简化系统的数学表达式。一般选择各个微分方程的输出及其各阶导数，也称相变量，作为状态变量，利用矩阵形式来表达系统，本质上依然是一系列微分方程。在状态空间模型基础上利用相应的方法可以对这些微分方程求解以获得汽轮机调节系统对输入量的动态响应。

1.1. 汽轮机转速调节系统运动的微分方程
1.1.1 汽轮发电机转子的运动方程
汽轮发电机组运行时，作用在转子上的力矩有三个：一是汽轮机的蒸汽主动力矩M_T，二是发电机反动力矩M_G，三是摩擦力矩M_F。由于摩擦力矩与蒸汽主动力矩和发电机反动力矩相比非常小，常常可以忽略不计。当M_T和M_G有变化时，转子的转速就会发生变化。若令蒸汽主动力矩增量为image，发电机反动力矩增量为image，相应的转动角速度增量为Δw，J为转动惯量，可得以下数学表达式：
image(1-1)
当在额定工况下时，额定蒸汽力矩等于额定反动力矩，即M_T0=M_G0。下标“0”表示额定值。将式(1-1)两端同时除以额定蒸汽力矩M_T0可得：
image(1-2)
考虑流进通流部分的蒸汽量与蒸汽汽室中的压力成正比，整理后得：
image(1-3)
式中：image为转子飞升时间常数；image为转速变化相对值；image为蒸汽汽室中压力变化相对值；β为汽轮发电机组自平衡系数；ψ为外界负荷扰动变化的相对值，其方向规定为增加负荷为正扰动。当自平衡作用很小时β-->0，式(1-3)可改写为：
image(1-4)
1.1.2 蒸汽容积方程
当进入汽轮机蒸汽汽室的蒸汽流量D1和流出的蒸汽流量D2有变化时就会引起汽室中蒸汽比重的变化，可表示为：
image(1-5)
式中：V为蒸汽汽室及管道的总容积；γ为蒸汽比重。这个容积方程和转子方程的形式相似，同理经转换可得：
image(1-6)
式中：image为汽室中蒸汽压力变化的相对值；image、image为阀门开度变化的相对值；image为汽室及管路总容积中的蒸汽流量与额定流量之比，其物理意义是以额定流量D0入汽室V中，使比重从0升高到γ0所需要的时间为Tρ。对于汽轮机调节阀后的蒸汽喷嘴室，由于其出口为喷嘴组，面积是固定不变的，所以Δh2=0、ξ2=0，于是式(1-6)可写为：
image(1-7)
1.1.3 调速器方程
对于理想的调速器（调速器的质量好、摩擦阻力小），当转速变化时，调速器的输出（如滑环位移）能迅速跟上，可认为在时间上没有滞后。这样相应转速输入就有相应位移的输出。
对应于δ的转速变化相对值，滑环位移变化为ΔXmax，于是可得到调速器方程为：
image(1-8)
式中：image为调速器滑环的位移变化相对值。
1.1.4 油动机方程
对于断流式双侧进油油动机，当转速升高时，调速器滑环将错油门活塞向上移动一个距离，打开油动机窗口使动力油进入油动机活塞上部，活塞下部则与排油管路相通而排油，油动机活塞在上下油压差的作用下关小调节气阀。
当忽视油动机位移时的惯性力与活塞上下油压作用的面积差值，并根据活塞上下进出油量相等的条件，油动机的微分方程为：
image(1-9)
式中：image为油动机活塞位移变化的相对值；image为错油门位移变化相对值；image为油动机时间常数，其物理意义是错油门开度保持在image时油动机活塞从满载位置关小到空载位置所需要的时间。
1.1.5 错油门运动方程
调速器、错油门和油动机构成反馈机构，当转速升高时，调速器滑环位移上升引起错油门抬高，打开油口使动力油进入油动机活塞上部，活塞下部从排油管路排油，此时油动机下降，关小汽阀。由于滑环、错油门活塞杆顶点和油动机活塞杆顶点构成杠杆，油动机下降使得以滑环为支点，错油门也下降关闭油口，这是负反馈的过程。所以在运动过程中错油门的位移即为滑环引起的位移和油动机引起的位移之差，可由下面的公式表示：
image(1-10)
1.2 汽轮机转速调节系统的状态空间模型
整理式(1-4)、(1-7)、(1-8)、(1-9)、(1-10)，同时注意以下两点：（1）因为油动机活塞杆直接与调节阀相连所以油动机相对位移μ与调节阀相对开度变化image应该相等，所以补充方程：image；（2）必须使系统闭环负反馈，例如转速增加时，油动机应下移，蒸汽流量应减小，所以调速器方程应改写为image，这样才反映了系统的真实情况并且是满足负反馈闭环系统要求；（3）从错油门运动方程推导过程可知，调速系统动作结束后，滑环位置没有回到原来的地方，代表着不能保持汽轮机原来的转速不变。这是由于调速系统稳定时错油门活塞必须处于原中间位置，即不变，而由于经过了调节过程，调节阀的开度必定发生变化，也就是油动机活塞在不同位置，根据杠杆原理，调速器滑环一定在不同的位置，与此相对应的转速数值必定与原来不同。这种调节系统动作后，稳定转速不能维持不变的的调节是有差调节，为了保证无差，需要给错油门额外的控制输入，这个控制输入可以是3.2节中MPC计算出的值，也可以是3.3节中PID控制器中的积分环节和微分环节的值。所以错油门方程应改写成image。最终，整个调节系统的微分方程组可表示为：
image(1-11)
选取φ、ρ、μ为状态变量，可得到状态方程为：
image(1-12)
将式(1-12)写出矩阵形式：
image(1-13)
输出即为转速变化相对值，输出方程为：
image(1-14)

2 基于模型预测控制的汽轮机转速控制

2.1 最优控制问题
最优控制问题就是研究在一定的约束条件下，如何选择控制规律u(t)才能达到最优的系统表现，通常系统的表现是综合分析的结果，没有绝对的最优，只有具体问题具体分析，在具体条件下达成的某种意义下的最优。比如考虑一个单输入单输出的系统（SISO），输出为y，要求其输出能跟踪预设的参考值r，误差可以表示为e=y-r，那么最优控制的目标是：
image(2-1)
式(2-1)是从输出的角度描述最优，所得结果越小，输出跟踪设定值的效果越好。除此之外，还可以从输入u的角度考虑，得到最优控制目标是：
image(2-2)
式(2-2)代表希望输入量 越小越好，从物理意义上来说就是能耗最低。两者结合，既希望输出跟踪设定值同时希望输入能耗少，可以得到最优控制的目标是：
image(2-3)
式中q和p是可调权重参数，跟据控制的重点不同，可以用于调节输出跟踪效果和输入低能耗之间的比重关系。
对于多输入多输出的系统（MIMO）image，状态变量是x，输入有r个，输出有m个，提出对输出和输入都有要求的二次型目标函数：
image(2-4)
式中R为r阶正定对称矩阵，Q为m阶正定对称矩阵。突出了系统要求有限时间内动态误差小，消耗的控制能量少。
以上的所有目标函数都是二次型函数，而二次型函数可以通过二次规划（Quadratic Programming）来求解最优控制序列。

2.2 模型预测控制原理
模型预测控制（MPC，Model Predictive Control）的基本思想是通过建立一个系统的动态模型，并在每一个控制时刻使用这个模型来预测系统未来的行为。最优化预测的系统行为，可以生成一个优化控制序列，然后通过执行第一个控制动作来调整系统状态，接着在下一个时刻重新计算和执行。这个过程反复进行，以使系统能够在未来的一段时间内持续优化一个特定的性能指标。
离散的模型方便计算机递推求解，已知系统的离散状态空间模型为：
image(2-5)
在当前k时刻预测未来的状态，预测区间为N。已知k时刻的状态x(k|k)，括号中后面的k代表要在k时刻进行预测。那么预测未来N个时间的状态x(k+1|k) ~ x(k+N|k)只由输入u(k|k) ~ u(k+N-1|k)决定。所以MPC控制器就是要找到N个控制输入使得预测区间内的状态满足要求的最优性能，一般希望区间内的输出跟踪设定值的误差越小越好，同时消耗的能量即输入的大小越小越好，目标函数一般为：
image(2-6)
然后选择这N个控制动作中的第一个控制动作作为当前k时刻的控制动作。最后把k+1时刻作为当前k时刻，再进行相同的操作。
只使用第一个控制动作，而舍弃其他的控制动作，是因为实际系统不是线性定常系统，可能发生一些变化，而且建模得出的模型不一定完全准确。所以在做完第一个动作后系统并不会像预测的那样转换到预测出的下一个状态。所以需要以k+1时刻测得的实际状态作为k+1时刻的预测区间N内预测所需的状态来校正计算，这个过程称为滚动优化控制（Receding Horizon Control）。
综上，MPC包括以下四个基本步骤：
（1）建立描述系统动态行为的数学模型，通常是差分方程，方便计算机迭代计算。
（2）在当前时刻基于已知的系统状态和待决策的控制输入，使用模型预测未来一段时间内的系统响应。
（3）基于预测的系统响应，通过求解一个与未来一段时间内的状态和控制量都有关的优化问题来计算最优的控制输入序列，以最大化或最小化一个性能指标（如累计误差、能耗等）。
（4）执行优化得到的控制输入序列中的第一个值，并将下一时刻的系统状态给下一个控制周期。
（5）返回第（2）步循环往复。

2.3 基于汽轮机转速控制的MPC控制器实现
为了能求解MPC最优控制规律，需要将其最优化目标函数转换成二次规划（Quadratic Programming）的形式，利用matlab中对于二次规划已经成熟的求解器来求解。二次规划的一般形式为：
image(2-7)
式中H是正定对称矩阵，f和x是向量。包含线性的等式约束和不等式约束。
2.3.1 目标函数
考虑在阶跃扰动下，希望汽轮机转速能够在控制规律作用下尽量保持不变。因为在定义转速控制系统的状态变量时，使用的是状态变化的相对值，所以在初始稳态所有状态值都是0，而且扰动出现时，希望作为输出的那个状态也就是转速变化相对值保持0不变，即设定值为0。每次计算最优问题时，考虑预测区间内的最优化目标函数应包括希望输出尽快跟踪设定值、希望输入能耗少、希望区间的终端输出越接近设定值越好。改写式(2-6)得目标函数为：
image(2-8)
式中x(k+i|k)是n个状态变量组成的列向量，u(k+i|k)是r个输入变量组成的列向量，则image、image、image都是对角权重矩阵，对角线上的元素分别表示运动过程中对各个状态大小的权重，对输入大小的权重和对终端时刻各个状态大小的权重。
把式(2-8)展开为：
image(2-9)
把式(2-9)写成矩阵形式，最终目标函数表示为：
image(2-10)
式中：
image
2.3.2 二次规划型目标函数
把式(2-10)所示的目标函数应用到MPC控制的具体案例时，可以转换为一般二次规划的形式。
假设预测区间为N，控制输入个数为P1，干扰输入个数为P2。把式(1-13)所示的汽轮机转速控制系统的模型离散化，如式(2-11)所示：
image(2-11)
式中 和 为原连续系统的系统矩阵和输入矩阵离散化的结果。
已知控制输入只有一个，作为错油门的输入，所以P1=1，对应输入矩阵为B矩阵第一列；干扰输入只有一个，是外界负荷干扰，进入汽轮发电机转子运动方程，所以P2=1，对应输入矩阵为B矩阵第二列。P1和P2在本汽轮机转速调节案例中不变。
在k时刻，预测N个步长的系统状态为：
image(2-12)
写成矩阵形式为：
image(2-13)
将其写为：
image(2-14)

式中：
image，image，image
image，image，image
把式(2-14)代入式(2-10)，得到：
image(2-15)
目标函数J的值是一个标量，所以式(2-15)中互为转置的相同下划线的项可以合并并把决策变量Uk提出，剩余项均为常数，对目标函数的最优化无影响。最终式(2-15)的形式与式(2-7)类似，所以可以用二次规划来求解。
最后考虑约束条件，输入是给错油门的，类似于调速器的滑环位移变化相对值，变化范围应在[-1,1]之间，所以约束条件为：1≤Uk≤1。

3 算例仿真

3.1 参数设置与仿真程序
设置转子飞升时间TM=10s，汽室及管路总容积中的蒸汽流量与额定流量之比Tρ=0.2s，使滑环位移变化为ΔXmax的转速变化相对值δ=0.05，油动机时间常数Ts=0.2s，得到汽轮机调节系统的状态空间表达式为：
image(3-1)
设置离散化采样时间T=0.01s，预测区间N=10，外界负荷扰动变化的相对值ψ=0.05，即5%额定负荷阶跃扰动，此时允许转速变化为±15转/分钟，相对额定转速3000转/分钟为±0.005，然后用MATLAB编写代码进行MPC仿真，用SIMULINK模块实现PID控制作为对比。
3.1.1 MPC程序
编写的MPC程序有一个主程序MPC.m进行参数设置，调用自定义的函数和画图；一个MPC_Matrices.m函数用于计算式(2-15)最终结果所需的各种矩阵，以便调用二次规划求解函数；一个Prediction.m函数用于预测区间内系统状态，给出最优控制序列，本质上是调用二次规划求解函数，然后返回第一个控制动作。
MPC.m如附录一，MPC_Matrices.m如附录二，Prediction.m如附录三。
3.1.2 PID仿真模块
用SIMULINK实现PID控制，模块和参数设置如图3-1和3-2所示。式(1-11)中image是把转速偏差反馈给系统，是比例控制，这也从另一个方面解释了1.2节中的有差控制。此时式(1-11)中image中的u这个控制输入就是积分环节和微分环节的值，而比例环节的值ξ可以和u加一起作为控制输入，这样可以把PID独立出来，那么原被控对象的系统矩阵中第三行第四列的元素就变为0了。
image

图3-1　PID控制模块图

image

图3-2　PID控制中汽轮机转速调节系统被控对象参数设置图

3.2 MPC控制
3.2.1 输出权重和输入权重
只关注作为输出的转速变化相对值，也就是第一个状态和输入大小时，根据式(2-8)取权重矩阵image，image，image时控制效果最好，转速变化相对值最大偏差为-0.0002，在0.4s时达到无差，如图3-3所示：
image

图3-3　最佳控制效果
当减小F矩阵以减小对终端输出的权重时，设置image，控制效果变差，系统不断振荡，难以稳定，如图3-4所示：
image

图3-4　减小F矩阵时的控制效果
在减小F矩阵情况下增加Q矩阵以增加对区间内输出大小的权重，设置image，控制效果回好，与图3-3的效果类似，转速变化相对值最大偏差为-0.0002，在0.6s时达到无差，如图3-5所示：
image

图3-5　减小F矩阵然后增加Q矩阵时的控制效果
在控制效果最佳的权重矩阵image，image时增加R，取R=10，更看重输入的能量大小，仿真结果显示输入值减小，但控制效果也变差，如图3-6所示：
image

图3-6　最佳Q和F矩阵下R=10时的控制效果

3.2.2 其他状态权重
设置image，image，R=0.1，即3.1.1节中只关注输出状态时的最佳的权重矩阵，得到转速变化相对值φ、汽室中蒸汽压力变化的相对值ρ、油动机活塞位移变化的相对值μ三个状态的变化如图3-7所示：
image

图3-7　最佳权重下各状态变化情况
设置image，image，此时比较看重ρ这个状态的大小，仿真结果显示ρ值明显减小，但是对转速的控制效果也变差，最大动态偏差变大，无法达到无差。三个状态的相应变化如图3-8所示：
image

图3-8　增加对蒸汽压力相对值的权重时各状态变化情况
从式(1-4)可以看出，在负荷扰动相对值ψ阶跃变化的情况下，要想汽轮机转速在稳态时保持不变，蒸汽压力变化的相对值ρ必须要抵消扰动的影响，那么就必定会有一定的超调使导数为正，使转速回升，最终🥩的大小与扰动值相等。而在目标函数的设计过程中的想要最小化的偏差全部都是与0相比的偏差，如果过分关注ρ和μ的大小，使其超调不足就一定会影响到对转速输出的控制。
设置image，image，此时比较看重油动机活塞位移变化的相对值μ的大小，仿真结果显示μ值明显减小，但同样对转速的控制效果也变差，最大动态偏差变大，同样无法达到无差。三个状态的相应变化如图3-9所示：
image

图3-9　增加对油动机活塞位移相对值的权重时各状态变化情况

3.3 PID控制与对比
设置PID控制器中Kp=20，Ki=10，Kd=20，仿真结果显示转速变化相对值最大偏差为-0.00008，相比MPC控制变化更小，但是调节速度缓慢，需要20s达到无差，如图3-10所示：
image

图3-10　PID控制
把MPC控制与PID控制的转速变化进行对比，虽然MPC控制的动态偏差更大，但是依然在允许范围内，PID控制的更加平滑，动态偏差小，但是响应速度非常缓慢。如图3-11所示：
image

图3-11　MPC与PID控制对比

3.4 模型不确定性下的MPC性能
式(1-4)是假设汽轮发电机组的自平衡作用很小，ρ-->0时转子的运动方程。如果实际的系统自平衡作用不可忽略或者自平衡性能有所改变时，MPC利用滚动优化控制，把k+1时刻测得的实际状态作为这个时刻预测所需的初始状态来校正计算，能承受多大的汽轮发电机组自平衡大小方面的模型不确定性。
考虑式(1-3)所示的含有汽轮发电机组自平衡作用的运动方程，可以得出实际系统的状态空间模型与机理建模的模型不同之处在于系统矩阵的第一个元素，应为image，即image。实际系统的状态空间模型如下：
image(3-2)
在最佳权重矩阵image，image，R=0.1设置β=10，MPC依然可以控制的很好，仿真结果如图3-12所示：
image

图3-12　自平衡系数β=10
设置β=19，转速变化相对值剧烈抖动，但是最后很快达到无差，如图3-13所示：
image

图3-13　自平衡系数β=19
当β=20，转速变化相对值开始等幅振荡，MPC无法控制，如图3-14所示。此时增加image，增加image，转速输出的振荡幅度有所下降，但频率依然很高，而且相比理想的控制效果相差很远，如图3-15所示。
image

图3-14　自平衡系数β=20
image

图3-15　自平衡系数β=20，增加权重矩阵Q和F
所以在汽轮发电机组的自平衡作用不确定下，MPC利用滚动优化控制可以适应β<20的不确定度。当不确定度很大时再增加对输出偏差的权重，起到的效果有限，因为二次规划型目标函数使用的依然是机理法建立的很不准确的模型，在计算出的控制输入序列作用下，实际状态与预测状态相差很远。

4 结论
本文基于模型预测控制（MPC）设计了汽轮机转速调节器，首先通过机理法建立汽轮机转速调节系统的状态空间模型，然后验证预测控制可以采用二次规划来求解最优控制序列，最后设置模型参数进行仿真。仿真结果表明，MPC控制器在权重参数合理设置下能够有效抑制负荷扰动引起的转速波动，动态响应迅速，动态偏差在可接受范围内。与PID控制相比，MPC具有更快的调节速度，虽然动态偏差略大，但仍满足允许范围，而PID控制的调节时间过于缓慢。综合而言，MPC相比PID控制有一定优势，为汽轮机转速调节提供了一种有效的优化方法。然而，本文设计的MPC也有一定局限性：
（1）依赖精确的模型，如果模型非常不精确则预测出的状态和求得的控制输入序列就缺乏参考意义。
（2）难以实时控制，因为完成k时刻的最优化计算需要当前的状态量，然后经过算法解算得到最优控制序列，导致控制不连续不实时。
（3）如果要对设定值阶跃的响应进行控制，需要重新设计目标函数，因为此时输出跟踪的目标不再是0。
（4）需要干扰信号已知或者能实时测量干扰。

附录
附录1：MPC.m程序

点击查看代码 
clear ; close all; clc;
%% 第一步，定义状态空间矩阵
A = [0 0.1 0; 
     0 -5  5;
     -100  0 -5];
B = [0 -0.1; 
     0  0;
     5  0];
C = [1 0 0];
sysd = c2d(ss(A,B,C,0),0.01);%T=0.01s
A = sysd.A;
B = sysd.B;

n= size (A,1);
p = 1;%作为控制的输入只有一个

%% 定义Q矩阵，n x n 矩阵
Q=[10 0 0;
    0  0 0;
    0  0 0];

%% 定义F矩阵，n x n 矩阵
F=[100000000 0 0;
    0  0 0;
    0  0 0];

%% 定义R矩阵，p x p 矩阵
R=0.1;

%% 定义step数量k
k_steps=150;    %仿真时间k_steps*0.01s

%% 定义矩阵 X_K， n x k 矩 阵
X_K = zeros(n,k_steps);

%% 初始状态变量值， n x 1 向量
X_K(:,1) =[0;
           0;
           0];

%% 定义要决策的输入矩阵 U_K，即只有一个控制输入(p=1)，另一个是干扰。 p x k 矩阵
U_K=[zeros(p,k_steps);
      0.05*ones(1,k_steps)];% 5%额定负荷阶跃扰动

%% 定义预测区间K
N = 10;

%% Call MPC_Matrices 函数 求得 E,H矩阵 
[E,H,L]=MPC_Matrices(A,B,Q,R,F,N,p);

%% 计算每一步的状态变量的值
for k = 1 : k_steps 
%% 求得U_K(:,k)
U_K(1,k) = Prediction(X_K(:,k),E,H,L,N,p);
%% 计算第k+1步时状态变量的值
X_K(:,k+1) = (A+[ 0 0 0;0 0 0; 0 0 0])*X_K(:,k) + B(:,1)*U_K(1,k) + B(:,2)*U_K(2,k);
end

%% 绘制状态变量和输入的变化
subplot(2, 1, 1);
hold;
for i =1 %:size (X_K,1)
plot (X_K(i,:));
end
legend('\phi','\rho','\mu')
xlabel('时间（*0.01s）')
ylabel('变化相对值')
hold off;

%%输入的图

subplot (2, 1, 2);
hold;
for i =1 : size (U_K,1)
plot (U_K(i,:));
xlabel('时间（*0.01s）')
ylabel('变化相对值')
end
legend('u','\psi')

%PID对比

% hold;
% for i =1 %:size (X_K,1)
% plot (X_K(i,:));
% end
% plot(100*PID.time,PID.signals.values)
% xlim([-50,2500])
% legend('MPC','PID')
% xlabel('时间（*0.01s）')
% ylabel('转速变化相对值')
% hold off;

%% x1 x2 x3分开画
% subplot(2, 1, 1);
% hold;
% for i =1
% plot (X_K(i,:));
% end
% legend('\phi')
% xlabel('时间（*0.01s）')
% ylabel('变化相对值')
% hold off;
% 
% subplot(2, 1, 2);
% hold;
% for i =2 :size (X_K,1)
% plot (X_K(i,:));
% end
% legend('\rho','\mu')
% xlabel('时间（*0.01s）')
% ylabel('变化相对值')
% hold off;

附录2：MPC_Matrices.m程序

点击查看代码 
function  [E , H,L]=MPC_Matrices(A,B,Q,R,F,N,p)
n=size(A,1);   % A 是 n x n 矩阵, 得到 n
% p=size(B,2);   % B 是 n x p 矩阵, 得到 p

% 定义M 和 C 
M=[eye(n);zeros(N*n,n)]; % 初始化 M 矩阵. M 矩阵是 (N+1)n x n的， 
                         % 它上面是 n x n 单位阵;, 这一步先把下半部
                         % 分写成 0 
C=zeros((N+1)*n,N*p); % 初始化 C 矩阵, 这一步令它有 (N+1)n x NP 个 0
D=zeros((N+1)*n,N*p);
pussy=0.05*ones(N*p,1);
tmp=eye(n);  %定义一个n x n 的 I 矩阵

%　更新Ｍ和C
for i=1:N % 循环，i 从 1到 N
    rows =i*n+(1:n); %定义当前行数，从i x n开始，共n行 
    C(rows,:)=[tmp*B(:,1),C(rows-n, 1:end-p)]; %将c矩阵填满
    D(rows,:)=[tmp*B(:,2),D(rows-n, 1:end-p)]; %将d矩阵填满
    tmp= A*tmp; %每一次将tmp左乘一次A
    M(rows,:)=tmp; %将M矩阵写满
end 

% 定义Q_bar和R_bar
Q_bar = kron(eye(N),Q);
Q_bar = blkdiag(Q_bar,F);
R_bar = kron(eye(N),R); 

% 计算G, E, H
G=M'*Q_bar*M; % G: n x n
E=M'*Q_bar*C; % E: NP x n
L=pussy'*D'*Q_bar*C;
H=C'*Q_bar*C+R_bar; % NP x NP 
end

附录3：Prediction.m程序

点击查看代码 
function u_k= Prediction(x_k,E,H,L,N,p)
U_k = zeros(N*p,1); % NP x 1
U_k = quadprog(2*H,2*(x_k'*E+L)',[],[],[],[],-1*ones(N,1),ones(N,1));
u_k = U_k(1,1); % 取第一个结果
end
posted @ 2026-03-14 16:34  -86CD  阅读(0)  评论(1)    收藏  举报