永磁同步电机控制系统——模型预测控制（MPC）

1 MPC原理

模型预测控制（Model Predictive Control, MPC）是近年来被广泛讨论的反馈控制策略。模型预测控制的机理可描述为：在每一采样时刻，根据获得的当前测量信息，在线求解一个有限时域开环优化问题，并将得到的控制序列的第一个元素作用于被控对象。在下一个采样时刻，重复上述过程，用新的测量值刷新优化问题并重新求解^1。

MPC与传统控制方法相比的几大优势

MPC在线求解开环优化问题，获得最优控制律。而传统控制方法通常离线求解一个反馈控制律，并将其一直作用于系统。MPC在处理非线性、等式约束和不等式约束问题上具有较大优势。

2 MPC求解无约束问题步骤

为引入积分以减小或者消除静态误差，使用增量式状态空间模型：

Δx(k+1)=AΔx(k)+BuΔu(k)+BdΔd(k)(1)\Delta x(k+1)=A\Delta x(k)+B_u\Delta u(k)+B_d\Delta d(k) \tag {1} Δx(k+1)=AΔx(k)+Bu​Δu(k)+Bd​Δd(k)(1)

yc(k)=CcΔx(k)+yc(k−1)(2)y_c(k)=C_c\Delta x(k)+y_c(k-1) \tag 2 yc​(k)=Cc​Δx(k)+yc​(k−1)(2)

在当前时刻kkk，测量值为x(k)x(k)x(k)，计算Δx(k)=x(k)−x(k−1)\Delta x(k)=x(k)-x(k-1)Δx(k)=x(k)−x(k−1)，将其作为预测系统未来动态的起点，预测未来ppp步的状态增量。

Δx(k+1∣k)=AΔx(k)+BuΔu(k)+BdΔd(k)Δx(k+2∣k)=AΔx(k+1∣k)+BuΔu(k+1)+BdΔd(k+1)=A2Δx(k)+ABuΔu(k)+BuΔu(k+1)+ABdΔd(k)⋮Δx(k+p∣k)=AΔx(k+p−1∣k)+BuΔu(k+p−1)+BdΔd(k+p−1)=ApΔx(k)+Ap−1BuΔu(k)+Ap−2BuΔu(k+1)+...+Ap−mBuΔu(k+m−1)+Ap−1BdΔd(k)(3)\begin{aligned} \Delta x(k+1|k)=&A\Delta x(k)+B_u\Delta u(k)+B_d\Delta d(k)\\ \Delta x(k+2|k)=&A\Delta x(k+1|k)+B_u\Delta u(k+1)+B_d\Delta d(k+1)\\ =&A^2\Delta x(k)+AB_u\Delta u(k)+B_u\Delta u(k+1)+AB_d\Delta d(k)\\ &\vdots\\ \Delta x(k+p|k)=&A\Delta x(k+p-1|k)+B_u\Delta u(k+p-1)+B_d\Delta d(k+p-1)\\ =&A^p\Delta x(k)+A^{p-1}B_u\Delta u(k)+A^{p-2}B_u\Delta u(k+1)\\ &+...+A^{p-m}B_u\Delta u(k+m-1)+A^{p-1}B_d\Delta d(k) \end{aligned} \tag 3 Δx(k+1∣k)=Δx(k+2∣k)==Δx(k+p∣k)==​AΔx(k)+Bu​Δu(k)+Bd​Δd(k)AΔx(k+1∣k)+Bu​Δu(k+1)+Bd​Δd(k+1)A2Δx(k)+ABu​Δu(k)+Bu​Δu(k+1)+ABd​Δd(k)⋮AΔx(k+p−1∣k)+Bu​Δu(k+p−1)+Bd​Δd(k+p−1)ApΔx(k)+Ap−1Bu​Δu(k)+Ap−2Bu​Δu(k+1)+...+Ap−mBu​Δu(k+m−1)+Ap−1Bd​Δd(k)​(3)

进一步，由输出方程可以预测未来ppp步被控输出

yc(k+1∣k)=CcΔx(k+1∣k)+yc(k)=CcAΔx(k)+CcBuΔu(k)+CcBdΔd(k)+yc(k)⋮yc(k+p∣k)=CcΔx(k+p∣k)+yc(k+p−1∣k)=∑i=1pCcAiΔx(k)+∑i=1pCcAi−1BuΔu(k)+∑i=1p−1CcAi−1BuΔu(k+1)+...+∑i=1p−m+1CcAi−1BuΔu(k+m−1)+∑i=1pCcAi−1BdΔd(k)+yc(k)(4)\begin{aligned} y_c(k+1|k)=&C_c\Delta x(k+1|k)+y_c(k)\\ =&C_cA\Delta x(k)+C_cB_u\Delta u(k)+C_cB_d\Delta d(k)+y_c(k)\\ &\vdots\\ y_c(k+p|k)=&C_c\Delta x(k+p|k)+y_c(k+p-1|k)\\= &\sum_{i=1}^pC_cA^i\Delta x(k)+\sum_{i=1}^pC_cA^{i-1}B_u\Delta u(k)\\ &+\sum_{i=1}^{p-1}C_cA^{i-1}B_u\Delta u(k+1)+...\\ &+\sum_{i=1}^{p-m+1}C_cA^{i-1}B_u\Delta u(k+m-1)\\ &+\sum_{i=1}^{p}C_cA^{i-1}B_d\Delta d(k)+y_c(k) \end{aligned} \tag 4 yc​(k+1∣k)==yc​(k+p∣k)==​Cc​Δx(k+1∣k)+yc​(k)Cc​AΔx(k)+Cc​Bu​Δu(k)+Cc​Bd​Δd(k)+yc​(k)⋮Cc​Δx(k+p∣k)+yc​(k+p−1∣k)i=1∑p​Cc​AiΔx(k)+i=1∑p​Cc​Ai−1Bu​Δu(k)+i=1∑p−1​Cc​Ai−1Bu​Δu(k+1)+...+i=1∑p−m+1​Cc​Ai−1Bu​Δu(k+m−1)+i=1∑p​Cc​Ai−1Bd​Δd(k)+yc​(k)​(4)

定义ppp步预测输出向量和mmm步预测输入向量为：

Yp(k+1∣k)=[yc(k+1∣k)yc(k+2∣k)⋮yc(k+p∣k)]p×1,ΔU(k)=[Δu(k)Δu(k+1)⋮Δu(k+m−1)]m×1(5)Y_p(k+1|k)=\begin{bmatrix} y_c(k+1|k) \\ y_c(k+2|k)\\\vdots\\y_c(k+p|k) \\ \end{bmatrix}_{p\times1}, \Delta U(k)=\begin{bmatrix}\Delta u(k) \\ \Delta u(k+1)\\\vdots\\\Delta u(k+m-1) \\ \end{bmatrix}_{m\times1} \tag 5 Yp​(k+1∣k)=⎣⎢⎢⎢⎡​yc​(k+1∣k)yc​(k+2∣k)⋮yc​(k+p∣k)​⎦⎥⎥⎥⎤​p×1​,ΔU(k)=⎣⎢⎢⎢⎡​Δu(k)Δu(k+1)⋮Δu(k+m−1)​⎦⎥⎥⎥⎤​m×1​(5)

注：矩阵的下标表示矩阵中向量或标量的个数。例如，上式中，p×1p\times1p×1仅表示Yc(k+1∣k)Y_c(k+1|k)Yc​(k+1∣k)矩阵中ycy_cyc​的个数。

系统未来ppp步预测输出为：

Yp(k+1∣k)=SxΔx(k)+Iyc(k)+SdΔd(k)+SuΔU(k)(6)Y_p(k+1|k)=S_x\Delta x(k)+\Iota y_c(k)+S_d\Delta d(k)+S_u\Delta U(k) \tag 6 Yp​(k+1∣k)=Sx​Δx(k)+Iyc​(k)+Sd​Δd(k)+Su​ΔU(k)(6)

其中，

Sx=[CcA∑i=12CcAi⋮∑i=1pCcAi]p×1,I=[Inc×ncInc×nc⋮Inc×nc]p×1,Sd=[CcBd∑i=12CcAi−1Bd⋮∑i=1pCcAi−1Bd]p×1,Su=[CcBu00⋯0∑i=12CcAi−1BuCcBu0⋯0⋮⋮⋮⋱⋮∑i=1mCcAi−1Bu∑i=1m−1CcAi−1Bu⋯⋯CcBu⋮⋮⋮⋱⋮∑i=1pCcAi−1Bu∑i=1p−1CcAi−1Bu⋯⋯∑i=1p−m+1CcAi−1Bu]p×m(7)S_x=\begin{bmatrix}C_cA \\ \sum_{i=1}^2C_cA^i\\\vdots\\\sum_{i=1}^pC_cA^i\end{bmatrix}_{p\times1}, \Iota=\begin{bmatrix}I_{n_c\times n_c} \\ I_{n_c\times n_c}\\\vdots\\I_{n_c\times n_c}\end{bmatrix}_{p\times1}, S_d=\begin{bmatrix}C_cB_d \\ \sum_{i=1}^2C_cA^{i-1}B_d\\\vdots\\\sum_{i=1}^pC_cA^{i-1}B_d\end{bmatrix}_{p\times1},\\ S_u=\begin{bmatrix}C_cB_u & 0 & 0 &\cdots&0\\ \sum_{i=1}^2C_cA^{i-1}B_u & C_cB_u & 0 & \cdots & 0\\\vdots&\vdots&\vdots&\ddots&\vdots\\\sum_{i=1}^mC_cA^{i-1}B_u & \sum_{i=1}^{m-1}C_cA^{i-1}B_u&\cdots&\cdots&C_cB_u\\\vdots&\vdots&\vdots&\ddots&\vdots\\\sum_{i=1}^pC_cA^{i-1}B_u & \sum_{i=1}^{p-1}C_cA^{i-1}B_u&\cdots&\cdots&\sum_{i=1}^{p-m+1}C_cA^{i-1}B_u\end{bmatrix}_{p\times m} \tag {7} Sx​=⎣⎢⎢⎢⎡​Cc​A∑i=12​Cc​Ai⋮∑i=1p​Cc​Ai​⎦⎥⎥⎥⎤​p×1​,I=⎣⎢⎢⎢⎡​Inc​×nc​​Inc​×nc​​⋮Inc​×nc​​​⎦⎥⎥⎥⎤​p×1​,Sd​=⎣⎢⎢⎢⎡​Cc​Bd​∑i=12​Cc​Ai−1Bd​⋮∑i=1p​Cc​Ai−1Bd​​⎦⎥⎥⎥⎤​p×1​,Su​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​Cc​Bu​∑i=12​Cc​Ai−1Bu​⋮∑i=1m​Cc​Ai−1Bu​⋮∑i=1p​Cc​Ai−1Bu​​0Cc​Bu​⋮∑i=1m−1​Cc​Ai−1Bu​⋮∑i=1p−1​Cc​Ai−1Bu​​00⋮⋯⋮⋯​⋯⋯⋱⋯⋱⋯​00⋮Cc​Bu​⋮∑i=1p−m+1​Cc​Ai−1Bu​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​p×m​(7)

接下来，要明确优化问题的数学描述。我们希望被控输出接近参考输入，同时希望控制变化不要太大，建立如下目标函数：

J=∑i=1p∥Γy,iyc(k+i∣k)−r(k+1)∥2+∑i=1m∥Γu,iΔu(k+i−1)∥2(8)J=\sum_{i=1}^p\begin{Vmatrix}\Gamma_{y,i}y_c(k+i|k)-r(k+1)\end{Vmatrix}^2 +\sum_{i=1}^m\begin{Vmatrix}\Gamma_{u,i}\Delta u(k+i-1)\end{Vmatrix}^2 \tag 8 J=i=1∑p​∥∥​Γy,i​yc​(k+i∣k)−r(k+1)​∥∥​2+i=1∑m​∥∥​Γu,i​Δu(k+i−1)​∥∥​2(8)

其中，Γy,i=diag(Γy1,i,Γy2,i,⋯,Γync,i),Γu,i=diag(Γu1,i,Γu2,i,⋯,Γunu,i)\Gamma_{y,i}=diag(\Gamma_{y_{1},i},\Gamma_{y_{2},i},\cdots,\Gamma_{y_{nc},i}),\Gamma_{u,i}=diag(\Gamma_{u_{1},i},\Gamma_{u_{2},i},\cdots,\Gamma_{u_{nu},i})Γy,i​=diag(Γy1​,i​,Γy2​,i​,⋯,Γync​,i​),Γu,i​=diag(Γu1​,i​,Γu2​,i​,⋯,Γunu​,i​)。Γyj,i\Gamma_{y_{j},i}Γyj​,i​是在预测时刻iii对预测控制输出第jjj个分量误差的加权因子，加权因子越大，表明期望控制输出越接近给定的参考输入；Γuj,i\Gamma_{u_{j},i}Γuj​,i​是在预测时刻iii对控制增量第jjj个分量的加权因子，控制加权因子越大，表明期望的控制动作变化越小。

优化问题1：

min⁡ΔU(k)J(x(k),ΔU(k),m,p)(9)\min\limits_{\Delta U(k)}J(x(k),\Delta U(k),m,p) \tag 9 ΔU(k)min​J(x(k),ΔU(k),m,p)(9)

满足动力学方程(i=0,1,⋯,p)(i=0,1,\cdots,p)(i=0,1,⋯,p)：

Δx(k+i+1∣k)=AΔx(k+i∣k)+BuΔu(k+i)+BdΔd(k+i),Δx(k∣k)=Δx(k),yc(k+i∣k)=CcΔx(k+i∣k)+yc(k+i−1∣k)(10)\begin{aligned} \Delta x(k+i+1|k)&=A\Delta x(k+i|k)+B_u\Delta u(k+i)+B_d\Delta d(k+i),\\ \Delta x(k|k) &= \Delta x(k),\\ y_c(k+i|k)&=C_c\Delta x(k+i|k)+y_c(k+i-1|k) \end{aligned} \tag {10} Δx(k+i+1∣k)Δx(k∣k)yc​(k+i∣k)​=AΔx(k+i∣k)+Bu​Δu(k+i)+Bd​Δd(k+i),=Δx(k),=Cc​Δx(k+i∣k)+yc​(k+i−1∣k)​(10)

其中，

J(x(k),ΔU(k),m,p)=∥ΓyYp(k+1∣k)−R(k+1)∥2+∥ΓuΔU(k)∥2(11)\begin{aligned} J(x(k),\Delta U(k),m,p)=&\begin{Vmatrix}\Gamma_{y}Y_p(k+1|k)-R(k+1)\end{Vmatrix}^2 +\begin{Vmatrix}\Gamma_{u}\Delta U(k)\end{Vmatrix}^2 \end{aligned} \tag {11} J(x(k),ΔU(k),m,p)=​∥∥​Γy​Yp​(k+1∣k)−R(k+1)​∥∥​2+∥∥​Γu​ΔU(k)​∥∥​2​(11)

其中，Yp(k+1∣k)Y_p(k+1|k)Yp​(k+1∣k)由式(5)\href{#eq5}{(5)}(5)给出，加权矩阵为：

Γy=diag(Γy,1,Γy,2,⋯,Γy,p),Γu=diag(Γu,1,Γu,2,⋯,Γu,m)(12)\Gamma_{y}=diag(\Gamma_{y,1},\Gamma_{y,2},\cdots,\Gamma_{y,p}),\\ \Gamma_{u}=diag(\Gamma_{u,1},\Gamma_{u,2},\cdots,\Gamma_{u,m}) \tag {12} Γy​=diag(Γy,1​,Γy,2​,⋯,Γy,p​),Γu​=diag(Γu,1​,Γu,2​,⋯,Γu,m​)(12)

参考输入序列为：

R(k+1)=[r(k+1)r(k+2)⋮r(k+p)](13)R(k+1)=\begin{bmatrix}r(k+1)\\r(k+2)\\\vdots\\r(k+p)\end{bmatrix} \tag {13} R(k+1)=⎣⎢⎢⎢⎡​r(k+1)r(k+2)⋮r(k+p)​⎦⎥⎥⎥⎤​(13)

求解1：

为方便求解，定义辅助变量：

ρ=[ΓyYp(k+1∣k)−R(k+1)ΓuΔU(k)](14)\rho = \begin{bmatrix}\Gamma_{y}Y_p(k+1|k)-R(k+1)\\ \Gamma_{u}\Delta U(k) \end{bmatrix} \tag{14} ρ=[Γy​Yp​(k+1∣k)−R(k+1)Γu​ΔU(k)​](14)

则目标函数(11)\href{#eq11}{(11)}(11)改写为：

J(x(k),ΔU(k),m,p)=ρTρ(15)J(x(k),\Delta U(k),m,p)=\rho^T\rho \tag {15} J(x(k),ΔU(k),m,p)=ρTρ(15)

将预测方程(6)\href{#eq6}{(6)}(6)代入(14)\href{#eq14}{(14)}(14)：

ρ=[Γy(SxΔx(k)+Iyc(k)+SdΔd(k)+SuΔU(k))−R(k+1)ΓuΔU(k)]=[ΓySuΓu]ΔU(k)−[Γy(R(k+1)−SxΔx(k)−Iyc(k)−SdΔd(k)⏞Ep(k+1∣k))0]=[ΓySuΓu]ΔU(k)⏟Az−[Ep(k+1∣k)0]⏟b=Az−b(16)\begin{aligned} \rho = & \begin{bmatrix} \Gamma_{y}(S_x\Delta x(k)+\Iota y_c(k)+S_d\Delta d(k)+S_u\Delta U(k))-R(k+1)\\ \Gamma_{u}\Delta U(k) \end{bmatrix}\\ =&\begin{bmatrix}\Gamma_{y}S_u\\\Gamma_{u}\end{bmatrix}\Delta U(k)-\begin{bmatrix}\Gamma_{y}(\overbrace{R(k+1)-S_x\Delta x(k)-\Iota y_c(k)-S_d\Delta d(k)}^{E_p(k+1|k)})\\0\end{bmatrix}\\ =&\underbrace{\begin{bmatrix}\Gamma_{y}S_u\\\Gamma_{u}\end{bmatrix}\Delta U(k)}_{Az}-\underbrace{\begin{bmatrix}E_p(k+1|k)\\0\end{bmatrix}}_b\\ =&Az-b \end{aligned} \tag {16} ρ====​[Γy​(Sx​Δx(k)+Iyc​(k)+Sd​Δd(k)+Su​ΔU(k))−R(k+1)Γu​ΔU(k)​][Γy​Su​Γu​​]ΔU(k)−⎣⎡​Γy​(R(k+1)−Sx​Δx(k)−Iyc​(k)−Sd​Δd(k)​Ep​(k+1∣k)​)0​⎦⎤​Az[Γy​Su​Γu​​]ΔU(k)​​−b[Ep​(k+1∣k)0​]​​Az−b​(16)

其中，

z=ΔU(k)，A=[ΓySuΓu],b=[Ep(k+1∣k)0](17)z=\Delta U(k)， A=\begin{bmatrix}\Gamma_{y}S_u\\\Gamma_{u}\end{bmatrix}, b=\begin{bmatrix}E_p(k+1|k)\\0\end{bmatrix} \tag {17} z=ΔU(k)，A=[Γy​Su​Γu​​],b=[Ep​(k+1∣k)0​](17)

Ep(k+1∣k)=R(k+1)−SxΔx(k)−Iyc(k)−SdΔd(k)(18)E_p(k+1|k)=R(k+1)-S_x\Delta x(k)-\Iota y_c(k)-S_d\Delta d(k) \tag {18} Ep​(k+1∣k)=R(k+1)−Sx​Δx(k)−Iyc​(k)−Sd​Δd(k)(18)

因此优化问题1变为：

min⁡zρTρ,其中ρ=Az−b(19)\min\limits_{z}\rho^T\rho,其中\rho=Az-b \tag {19} zmin​ρTρ,其中ρ=Az−b(19)

极值条件为：

dρTρdz=2(dρdz)T,ρ=2AT(Az−b)=0(20)\frac{{\rm d}\rho^T\rho}{{\rm d}z}=2(\frac{{\rm d}\rho}{{\rm d}z})^T,\rho=2A^T(Az-b)=0 \tag {20} dzdρTρ​=2(dzdρ​)T,ρ=2AT(Az−b)=0(20)

得极值解为

z∗=(ATA)−1ATb(21)z^*=(A^TA)^{-1}A^Tb \tag {21} z∗=(ATA)−1ATb(21)

又由于

d2(ρTρ)dz2=2ATA>0(22)\frac{{\rm d}^2(\rho^T\rho)}{{\rm d}z^2}=2A^TA>0 \tag {22} dz2d2(ρTρ)​=2ATA>0(22)

则式(21)\href{#eq21}{(21)}(21)是取得最小值的解。将式(17)\href{#eq17}{(17)}(17)代入式(21)\href{#eq21}{(21)}(21)，得到优化问题1的最优解，即kkk时刻的最优控制序列为

ΔU∗(k)=(SuTΓyTΓySu+ΓuTΓu)−1SuTΓyTΓyEp(k+1∣k)(23)\Delta U^*(k)=(S_u^T\Gamma_y^T\Gamma_yS_u+\Gamma_u^T\Gamma_u)^{-1}S_u^T\Gamma_y^T\Gamma_yE_p(k+1|k) \tag {23} ΔU∗(k)=(SuT​ΓyT​Γy​Su​+ΓuT​Γu​)−1SuT​ΓyT​Γy​Ep​(k+1∣k)(23)

其中Ep(k+1∣k)E_p(k+1|k)Ep​(k+1∣k)由式(18)\href{#eq18}{(18)}(18)计算。

将最优控制序列的第一个元素将作用于系统，即

Δu(k)=[Inc×nc0⋯0]1×mΔU∗(k)=[Inc×nc0⋯0]1×m(SuTΓyTΓySu+ΓuTΓu)−1SuTΓyTΓyEp(k+1∣k)(24)\begin{aligned} \Delta u(k) =&\begin{bmatrix}I_{n_c\times n_c}& 0 & \cdots & 0\end{bmatrix}_{1\times m}\Delta U*(k)\\ =&\begin{bmatrix}I_{n_c\times n_c}& 0 & \cdots & 0\end{bmatrix}_{1\times m}(S_u^T\Gamma_y^T\Gamma_yS_u+\Gamma_u^T\Gamma_u)^{-1}S_u^T\Gamma_y^T\Gamma_yE_p(k+1|k) \end{aligned} \tag {24} Δu(k)==​[Inc​×nc​​​0​⋯​0​]1×m​ΔU∗(k)[Inc​×nc​​​0​⋯​0​]1×m​(SuT​ΓyT​Γy​Su​+ΓuT​Γu​)−1SuT​ΓyT​Γy​Ep​(k+1∣k)​(24)

定义预测控制增益

Kmpc=[Inc×nc0⋯0]1×m(SuTΓyTΓySu+ΓuTΓu)−1SuTΓyTΓy(25)K_{mpc}=\begin{bmatrix}I_{n_c\times n_c}& 0 & \cdots & 0\end{bmatrix}_{1\times m}(S_u^T\Gamma_y^T\Gamma_yS_u+\Gamma_u^T\Gamma_u)^{-1}S_u^T\Gamma_y^T\Gamma_y \tag {25} Kmpc​=[Inc​×nc​​​0​⋯​0​]1×m​(SuT​ΓyT​Γy​Su​+ΓuT​Γu​)−1SuT​ΓyT​Γy​(25)

那么控制增量可由下式计算

Δu(k)=KmpcEp(k+1∣k)(26)\Delta u(k)=K_{mpc}E_p(k+1|k) \tag {26} Δu(k)=Kmpc​Ep​(k+1∣k)(26)

算法总结：

初始化：设定预测时域p和控制时域m，初始值u(−1)=0u(-1)=0u(−1)=0，x(−1)=0x(-1)=0x(−1)=0，由式(7)\href{#eq7}{(7)}(7)计算SxS_xSx​，III，SdS_dSd​和SuS_uSu​，由式(25)\href{#eq25}{(25)}(25)计算KmpcK_{mpc}Kmpc​。k≥0k≥0k≥0时刻，得到测量值x(k)x(k)x(k)和d(k)d(k)d(k)。由式(2)\href{#eq2}{(2)}(2)计算yc(k)y_c(k)yc​(k)，计算Δx(k)=x(k)−x(k−1)\Delta x(k)=x(k)-x(k-1)Δx(k)=x(k)−x(k−1)。由式(18)\href{#eq18}{(18)}(18)计算误差 Ep(k+1∣k)E_p(k+1|k)Ep​(k+1∣k)。由式(26)\href{#eq26}{(26)}(26)计算控制变量增量 Δu(k)\Delta u(k)Δu(k)。将u(k)=u(k−1)+Δu(k)u(k)=u(k-1)+\Delta u(k)u(k)=u(k−1)+Δu(k)作用于系统。在k+1k+1k+1时刻，测量x(k+1)x(k+1)x(k+1)和d(k+1)d(k+1)d(k+1)的值，并且令k=k+1k=k+1k=k+1，返回第2步。 注：如果采用时不变加权因子，则KmpcK_{mpc}Kmpc​可离线计算，否则KmpcK_{mpc}Kmpc​需在循环内进行，即在线计算

3 基于MPC的电流环设计

永磁同步电机电压方程如下：

diddt=−RsLdid+ωeLqLdiq+udLddiqdt=−ωeLdLqid−RsLqiq+uqLq−ωeψmLq(27)\begin{aligned} \frac {{\rm d}i_d}{{\rm d}t}=& -\frac {R_s} {L_d}i_d+\frac {\omega_eL_q}{L_d}i_q+\frac {u_d}{L_d}\\ \frac {{\rm d}i_q}{{\rm d}t}=& -\frac {\omega_eL_d}{L_q}i_d -\frac {R_s} {L_q}i_q+\frac {u_q}{L_q}-\frac {\omega_e\psi_m}{L_q} \end{aligned} \tag {27} dtdid​​=dtdiq​​=​−Ld​Rs​​id​+Ld​ωe​Lq​​iq​+Ld​ud​​−Lq​ωe​Ld​​id​−Lq​Rs​​iq​+Lq​uq​​−Lq​ωe​ψm​​​(27)

离散化并改写成增量式

Δx(k+1)=AΔx(k)+BuΔu(k)(28)\Delta x(k+1)=A\Delta x(k)+B_u\Delta u(k) \tag {28} Δx(k+1)=AΔx(k)+Bu​Δu(k)(28)

其中，

x=[idiq],u=[uduq]A=[1−RsTsLdωeLqTsLd−ωeLdTsLq1−RsTsLq],Bu=[TsLdTsLq](29)x=\begin{bmatrix}i_d & i_q\end{bmatrix},u=\begin{bmatrix}u_d & u_q\end{bmatrix}\\ A=\begin{bmatrix}1-\frac {R_sT_s}{L_d} & \frac {\omega_eL_qT_s}{L_d} \\ -\frac {\omega_eL_dT_s}{L_q} & 1-\frac {R_sT_s}{L_q}\end{bmatrix}, B_u=\begin{bmatrix}\frac {T_s}{L_d} & \frac {T_s}{L_q}\end{bmatrix} \tag {29} x=[id​​iq​​],u=[ud​​uq​​]A=[1−Ld​Rs​Ts​​−Lq​ωe​Ld​Ts​​​Ld​ωe​Lq​Ts​​1−Lq​Rs​Ts​​​],Bu​=[Ld​Ts​​​Lq​Ts​​​](29)

控制输出

yc(k)=CΔx(k)+yc(k−1),C=[1001](30)y_c(k)=C\Delta x(k)+y_c(k-1),C=\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix} \tag {30} yc​(k)=CΔx(k)+yc​(k−1),C=[10​01​](30)

d轴电流给定为0，给定参考速度为100rad/s，负载转矩0.25Nm，matlab仿真结果如下。稳定时，d轴电流为0，q轴电流为7.3369A，机械角速度为100rad/s。

4 基于MPC的电流环和速度环设计

结合永磁同步电机电气方程和机械方程：

diddt=−RsLdid+pnωmLqLdiq+udLddiqdt=−pnωmLdLqid−RsLqiq+uqLq−pnωmψmLqdωmdt=3pnψm2Jiq+3pn2J(Ld−Lq)idiq−BJωm−TlJ(31)\begin{aligned} \frac {{\rm d}i_d}{{\rm d}t}=& -\frac {R_s} {L_d}i_d+\frac {p_n\omega_mL_q}{L_d}i_q+\frac {u_d}{L_d}\\ \frac {{\rm d}i_q}{{\rm d}t}=& -\frac {p_n\omega_mL_d}{L_q}i_d -\frac {R_s} {L_q}i_q+\frac {u_q}{L_q}-\frac {p_n\omega_m\psi_m}{L_q}\\ \frac{{\rm d}\omega_m}{{\rm d}t}=&\frac{3p_n\psi_m}{2J}i_q+\frac{3p_n}{2J}(L_d-L_q)i_di_q-\frac{B}{J}\omega_m-\frac{T_l}{J} \end{aligned} \tag {31} dtdid​​=dtdiq​​=dtdωm​​=​−Ld​Rs​​id​+Ld​pn​ωm​Lq​​iq​+Ld​ud​​−Lq​pn​ωm​Ld​​id​−Lq​Rs​​iq​+Lq​uq​​−Lq​pn​ωm​ψm​​2J3pn​ψm​​iq​+2J3pn​​(Ld​−Lq​)id​iq​−JB​ωm​−JTl​​​(31)

这里讨论表贴式永磁同步电机，即Ld=Lq=LsL_d=L_q=L_sLd​=Lq​=Ls​。

diddt=−RsLsid+pnωmiq+udLsdiqdt=−pnωmid−RsLsiq+uqLs−pnωmψmLsdωmdt=3pnψm2Jiq−BJωm−TlJ(32)\begin{aligned} \frac {{\rm d}i_d}{{\rm d}t}=& -\frac {R_s} {L_s}i_d+p_n\omega_mi_q+\frac {u_d}{L_s}\\ \frac {{\rm d}i_q}{{\rm d}t}=& -{p_n\omega_m}i_d -\frac {R_s} {L_s}i_q+\frac {u_q}{L_s}-\frac {p_n\omega_m\psi_m}{L_s}\\ \frac{{\rm d}\omega_m}{{\rm d}t}=&\frac{3p_n\psi_m}{2J}i_q-\frac{B}{J}\omega_m-\frac{T_l}{J} \end{aligned} \tag {32} dtdid​​=dtdiq​​=dtdωm​​=​−Ls​Rs​​id​+pn​ωm​iq​+Ls​ud​​−pn​ωm​id​−Ls​Rs​​iq​+Ls​uq​​−Ls​pn​ωm​ψm​​2J3pn​ψm​​iq​−JB​ωm​−JTl​​​(32)

令x=[idiqωm],u=[uduq]x=\begin{bmatrix}i_d & i_q & \omega_m\end{bmatrix},u=\begin{bmatrix}u_d & u_q \end{bmatrix}x=[id​​iq​​ωm​​],u=[ud​​uq​​]，则

x˙=[−RsLs000−RsLs−pnψmLs03pnψm2J−BJ]x+[−1Ls00−1Ls00]u+[pnωmiq−pnωmid0]⏟非线性项+[00−TlJ](33)\dot x=\begin{bmatrix}-\frac{R_s}{L_s} & 0 & 0\\0 & -\frac{R_s}{L_s}&-\frac{p_n\psi_m}{L_s}\\0&\frac{3p_n\psi_m}{2J}&-\frac{B}{J}\end{bmatrix}x +\begin{bmatrix}-\frac{1}{L_s} & 0 \\0 & -\frac{1}{L_s}\\0&0\end{bmatrix}u +\underbrace{\begin{bmatrix}p_n\omega_mi_q\\-p_n\omega_mi_d\\0\end{bmatrix}}_{非线性项} +\begin{bmatrix}0\\0\\-\frac{T_l}{J}\end{bmatrix} \tag {33} x˙=⎣⎢⎡​−Ls​Rs​​00​0−Ls​Rs​​2J3pn​ψm​​​0−Ls​pn​ψm​​−JB​​⎦⎥⎤​x+⎣⎡​−Ls​1​00​0−Ls​1​0​⎦⎤​u+非线性项⎣⎡​pn​ωm​iq​−pn​ωm​id​0​⎦⎤​​​+⎣⎡​00−JTl​​​⎦⎤​(33)

假设x˙=f(x,u)\dot x=f(x,u)x˙=f(x,u)，根据全微分方程Δx˙=∂f∂xΔx+∂f∂uΔu\Delta \dot x=\frac{\partial f}{\partial x}\Delta x+\frac{\partial f}{\partial u}\Delta uΔx˙=∂x∂f​Δx+∂u∂f​Δu，将上式线性化：

Δx˙=[−RsLspnωm,tpniq,t−pnωm,t−RsLs−pnψmLs−pnid,t03pnψm2J−BJ]Δx+[−1Ls00−1Ls00]Δu(34)\Delta \dot x=\begin{bmatrix}-\frac{R_s}{L_s} & p_n\omega_{m,t} & p_ni_{q,t}\\-p_n\omega_{m,t} & -\frac{R_s}{L_s}&-\frac{p_n\psi_m}{L_s}-p_ni_{d,t}\\0&\frac{3p_n\psi_m}{2J}&-\frac{B}{J}\end{bmatrix}\Delta x +\begin{bmatrix}-\frac{1}{L_s} & 0 \\0 & -\frac{1}{L_s}\\0&0\end{bmatrix}\Delta u \tag {34} Δx˙=⎣⎢⎡​−Ls​Rs​​−pn​ωm,t​0​pn​ωm,t​−Ls​Rs​​2J3pn​ψm​​​pn​iq,t​−Ls​pn​ψm​​−pn​id,t​−JB​​⎦⎥⎤​Δx+⎣⎡​−Ls​1​00​0−Ls​1​0​⎦⎤​Δu(34)

对式(34)\href{#eq34}{(34)}(34)离散化得到：

Δx(k+1)=AΔx(k)+BuΔu(k)(35)\Delta x(k+1)=A\Delta x(k)+B_u\Delta u(k) \tag {35} Δx(k+1)=AΔx(k)+Bu​Δu(k)(35)

其中

A=[1−RsLsTspnωm,tTspniq,tTs−pnωm,tTs1−RsLsTs−pnψmLsTs−pnid,tTs03pnψm2JTs1−BJTs],Bu=[−TsLs00−TsLs00](36)A=\begin{bmatrix}1-\frac{R_s}{L_s}T_s & p_n\omega_{m,t}T_s& p_ni_{q,t}T_s\\-p_n\omega_{m,t}T_s& 1-\frac{R_s}{L_s}T_s&-\frac{p_n\psi_m}{L_s}T_s-p_ni_{d,t}T_s\\0&\frac{3p_n\psi_m}{2J}T_s&1-\frac{B}{J}T_s\end{bmatrix}, B_u=\begin{bmatrix}-\frac{T_s}{L_s} & 0 \\0 & -\frac{T_s}{L_s}\\0&0\end{bmatrix} \tag {36} A=⎣⎢⎡​1−Ls​Rs​​Ts​−pn​ωm,t​Ts​0​pn​ωm,t​Ts​1−Ls​Rs​​Ts​2J3pn​ψm​​Ts​​pn​iq,t​Ts​−Ls​pn​ψm​​Ts​−pn​id,t​Ts​1−JB​Ts​​⎦⎥⎤​,Bu​=⎣⎡​−Ls​Ts​​00​0−Ls​Ts​​0​⎦⎤​(36)

控制输出

yc(k)=CΔx(k)+yc(k−1),C=[1001](37)y_c(k)=C\Delta x(k)+y_c(k-1),C=\begin{bmatrix}1 & 0 \\ 0 & 1\end{bmatrix} \tag {37} yc​(k)=CΔx(k)+yc​(k−1),C=[10​01​](37)

根据第2节中的算法步骤求解，得到控制变量

u(k)=Δu(k)+u(k−1)(38)u(k)=\Delta u(k)+u(k-1) \tag {38} u(k)=Δu(k)+u(k−1)(38)

d轴电流给定为0，给定参考速度为100rad/s，负载转矩0.25Nm，matlab仿真结果如下。稳定时，d轴电流为0，q轴电流为7.3366A，机械角速度为99.8453rad/s。

matlab仿真参考代码。

5 参考文献

陈虹. 模型预测控制[M]. 北京:科学出版社, . ↩︎