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)CcAΔx(k)+CcBuΔu(k)+CcBdΔd(k)+yc(k)⋮CcΔx(k+p∣k)+yc(k+p−1∣k)i=1∑pCcAiΔx(k)+i=1∑pCcAi−1BuΔu(k)+i=1∑p−1CcAi−1BuΔu(k+1)+...+i=1∑p−m+1CcAi−1BuΔu(k+m−1)+i=1∑pCcAi−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=⎣⎢⎢⎢⎡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=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡CcBu∑i=12CcAi−1Bu⋮∑i=1mCcAi−1Bu⋮∑i=1pCcAi−1Bu0CcBu⋮∑i=1m−1CcAi−1Bu⋮∑i=1p−1CcAi−1Bu00⋮⋯⋮⋯⋯⋯⋱⋯⋱⋯00⋮CcBu⋮∑i=1p−m+1CcAi−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,iyc(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)minJ(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)=∥∥ΓyYp(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} ρ=[ΓyYp(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)][ΓySuΓu]ΔU(k)−⎣⎡Γy(R(k+1)−SxΔx(k)−Iyc(k)−SdΔd(k)Ep(k+1∣k))0⎦⎤Az[ΓySuΓ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=[ΓySuΓ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变为:
minzρ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ΓySu+ΓuTΓu)−1SuTΓyTΓyEp(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×nc0⋯0]1×mΔU∗(k)[Inc×nc0⋯0]1×m(SuTΓyTΓySu+ΓuTΓu)−1SuTΓyTΓyEp(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×nc0⋯0]1×m(SuTΓyTΓySu+Γ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)=KmpcEp(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=−LdRsid+LdωeLqiq+Ldud−LqωeLdid−LqRsiq+Lquq−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=[idiq],u=[uduq]A=[1−LdRsTs−LqωeLdTsLdωeLqTs1−LqRsTs],Bu=[LdTsLqTs](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=[1001](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=−LdRsid+LdpnωmLqiq+Ldud−LqpnωmLdid−LqRsiq+Lquq−Lqpnωmψm2J3pnψmiq+2J3pn(Ld−Lq)idiq−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=−LsRsid+pnωmiq+Lsud−pnωmid−LsRsiq+Lsuq−Lspnωmψm2J3pnψmiq−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=[idiqωm],u=[uduq],则
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˙=⎣⎢⎡−LsRs000−LsRs2J3pnψm0−Lspnψm−JB⎦⎥⎤x+⎣⎡−Ls1000−Ls10⎦⎤u+非线性项⎣⎡pnωmiq−pnωmid0⎦⎤+⎣⎡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˙=⎣⎢⎡−LsRs−pnωm,t0pnωm,t−LsRs2J3pnψmpniq,t−Lspnψm−pnid,t−JB⎦⎥⎤Δx+⎣⎡−Ls1000−Ls10⎦⎤Δ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−LsRsTs−pnωm,tTs0pnωm,tTs1−LsRsTs2J3pnψmTspniq,tTs−LspnψmTs−pnid,tTs1−JBTs⎦⎥⎤,Bu=⎣⎡−LsTs000−LsTs0⎦⎤(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=[1001](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]. 北京:科学出版社, . ↩︎