1.一种分数阶模型预测控制的加热炉温度控制方法,其特征在于该方法的具体步骤是:
步骤1、建立实际过程中被控对象的扩展状态空间模型,具体是:
1.1采集实际过程对象的实时阶跃响应数据,利用该数据建立被控对象的分数阶传递函数模型,形式如下:
<mrow>
<mi>G</mi>
<mrow>
<mo>(</mo>
<mi>s</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<msup>
<mi>Ke</mi>
<mrow>
<mo>-</mo>
<mi>&tau;</mi>
<mi>s</mi>
</mrow>
</msup>
</mrow>
<mrow>
<msub>
<mi>c</mi>
<mn>1</mn>
</msub>
<msup>
<mi>s</mi>
<msub>
<mi>&alpha;</mi>
<mn>1</mn>
</msub>
</msup>
<mo>+</mo>
<msub>
<mi>c</mi>
<mn>0</mn>
</msub>
</mrow>
</mfrac>
</mrow>
其中,α1为微分阶次,c0,c1为相应的系数,s为拉普拉斯变换算子,K为模型增益,τ为模型的滞后时间;
1.2由Oustaloup近似方法得到微分算子sα的近似表达形式如下:
<mrow>
<msup>
<mi>s</mi>
<mi>&alpha;</mi>
</msup>
<mo>&ap;</mo>
<msub>
<mi>K</mi>
<mi>&alpha;</mi>
</msub>
<munderover>
<mo>&Pi;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mfrac>
<mrow>
<mi>s</mi>
<mo>+</mo>
<msubsup>
<mi>w</mi>
<mi>n</mi>
<mo>&prime;</mo>
</msubsup>
</mrow>
<mrow>
<mi>s</mi>
<mo>+</mo>
<msub>
<mi>w</mi>
<mi>n</mi>
</msub>
</mrow>
</mfrac>
</mrow>
其中,α为分数阶微分阶次,0<α<1,N为选定的近似阶次,wb和wh分别为选定的拟合频率的下限和上限;
1.3根据步骤1.2中的方法,将步骤1.1中的分数阶传递函数模型近似为整数阶高阶模型,进而将其在采样时间Ts下加零阶保持器离散化,得到如下形式的离散模型:
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>-</mo>
<msub>
<mi>F</mi>
<mn>1</mn>
</msub>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mi>F</mi>
<mn>2</mn>
</msub>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mn>...</mn>
<mo>-</mo>
<msub>
<mi>F</mi>
<msub>
<mi>L</mi>
<mi>S</mi>
</msub>
</msub>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<msub>
<mi>k</mi>
<mi>S</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<msub>
<mi>H</mi>
<mn>1</mn>
</msub>
<mi>u</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
<mo>-</mo>
<mi>d</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>H</mi>
<mn>2</mn>
</msub>
<mi>u</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>2</mn>
<mo>-</mo>
<mi>d</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mn>...</mn>
<mo>+</mo>
<msub>
<mi>H</mi>
<msub>
<mi>L</mi>
<mi>S</mi>
</msub>
</msub>
<mi>u</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<msub>
<mi>L</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<mi>d</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中,Fj,Hj(j=1,2,…,LS)均为离散近似后得到的系数,实际过程的时滞d=τ/Ts,LS为离散模型的长度,y(k)为k时刻的实际过程对象的模型输出,u(k-d-1)为实际过程对象在k-d-1时刻的输入值;
进一步将上述模型取一阶向后差分,得到如下形式:
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>&Delta;</mi>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>F</mi>
<mn>1</mn>
</msub>
<mi>&Delta;</mi>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>F</mi>
<mn>2</mn>
</msub>
<mi>&Delta;</mi>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mn>...</mn>
<mo>+</mo>
<msub>
<mi>H</mi>
<msub>
<mi>L</mi>
<mi>S</mi>
</msub>
</msub>
<mi>&Delta;</mi>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<msub>
<mi>L</mi>
<mi>S</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<msub>
<mi>H</mi>
<mn>1</mn>
</msub>
<mi>&Delta;</mi>
<mi>u</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
<mo>-</mo>
<mi>d</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>H</mi>
<mn>2</mn>
</msub>
<mi>&Delta;</mi>
<mi>u</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mn>2</mn>
<mo>-</mo>
<mi>d</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mn>...</mn>
<mo>+</mo>
<msub>
<mi>H</mi>
<msub>
<mi>L</mi>
<mi>S</mi>
</msub>
</msub>
<mi>&Delta;</mi>
<mi>u</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<msub>
<mi>L</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<mi>d</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中,Δ是差分算子;
1.4选取如下状态变量:
Δxm(k)=[Δy(k),Δy(k-1),…,Δy(k-LS+1),Δu(k-1),…,Δu(k-LS+1-d)]T
结合步骤1.3,得到被控对象的状态空间模型,形式如下:
Δxm(k+1)=AmΔxm(k)+BmΔu(k)
Δy(k+1)=CmΔxm(k+1)
其中,T为矩阵的转置符号,Δxm(k)的维数为(2LS+d-1)×1;
Bm=[0 … 0 1 0 … 0]T
Cm=[1 0 0 … 0 0 0 0]
1.5将步骤1.4中得到的状态空间模型转换成包含状态变量和输出跟踪误差的扩展状态空间模型,形式如下:
z(k+1)=Az(k)+BΔu(k)+CΔr(k+1)
其中,
<mrow>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>e</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&Delta;x</mi>
<mi>m</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>;</mo>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>e</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&Delta;x</mi>
<mi>m</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
e(k)=y(k)-r(k)
e(k+1)=e(k)+CmAmΔxm(k)+CmBmΔu(k)-Δr(k+1)
<mrow>
<mi>A</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mrow>
<msub>
<mi>C</mi>
<mi>m</mi>
</msub>
<msub>
<mi>A</mi>
<mi>m</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<msub>
<mi>A</mi>
<mi>m</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>;</mo>
<mi>B</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>C</mi>
<mi>m</mi>
</msub>
<msub>
<mi>B</mi>
<mi>m</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>B</mi>
<mi>m</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>;</mo>
<mi>C</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
r(k)为k时刻的跟踪设定值,e(k)为k时刻的输出误差,0为(2LS+d-1)×1维的零矩阵,A为(2LS+d)×(2LS+d)维矩阵,B,C均为(2LS+d)×1维矩阵;
步骤2、基于扩展状态空间模型设计被控对象的分数阶模型预测控制器,具体如下:
2.1预测未来k+i时刻模型输出的向量形式,
Z=Gz(k)+SΔU+ΨΔR
其中,
<mrow>
<mi>Z</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mi>P</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>;</mo>
<mi>G</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>A</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msup>
<mi>A</mi>
<mn>2</mn>
</msup>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msup>
<mi>A</mi>
<mi>P</mi>
</msup>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
ΔU=[Δu(k) Δu(k+1) … Δu(k+M-1)]T
ΔR=[Δr(k+1) Δr(k+2) … Δr(k+P)]T
r(k+i)=λiy(k)+(1-λi)c(k)
c(k)为k时刻的设定值,λ为柔化因子,P为预测时域,M为控制时域,r(k+i)为k+i时刻过程的跟踪设定值,i=1,2,…,P;
2.2选取被控对象的目标函数J,其形式如下:
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>J</mi>
<mo>=</mo>
<mmultiscripts>
<mi>I</mi>
<msub>
<mi>T</mi>
<mi>S</mi>
</msub>
<mrow>
<msub>
<mi>PT</mi>
<mi>S</mi>
</msub>
</mrow>
<msub>
<mi>&gamma;</mi>
<mn>1</mn>
</msub>
</mmultiscripts>
<mi>z</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mmultiscripts>
<mi>I</mi>
<msub>
<mi>T</mi>
<mi>S</mi>
</msub>
<mrow>
<msub>
<mi>MT</mi>
<mi>S</mi>
</msub>
</mrow>
<msub>
<mi>&gamma;</mi>
<mn>2</mn>
</msub>
</mmultiscripts>
<mi>&Delta;</mi>
<mi>u</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<msubsup>
<mo>&Integral;</mo>
<msub>
<mi>T</mi>
<mi>S</mi>
</msub>
<mrow>
<msub>
<mi>PT</mi>
<mi>S</mi>
</msub>
</mrow>
</msubsup>
<msup>
<mi>D</mi>
<mrow>
<mn>1</mn>
<mo>-</mo>
<msub>
<mi>&gamma;</mi>
<mn>1</mn>
</msub>
</mrow>
</msup>
<mi>z</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
<mi>z</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mi>d</mi>
<mi>t</mi>
<mo>+</mo>
<msubsup>
<mo>&Integral;</mo>
<msub>
<mi>T</mi>
<mi>S</mi>
</msub>
<mrow>
<msub>
<mi>MT</mi>
<mi>S</mi>
</msub>
</mrow>
</msubsup>
<msup>
<mi>D</mi>
<mrow>
<mn>1</mn>
<mo>-</mo>
<msub>
<mi>&gamma;</mi>
<mn>2</mn>
</msub>
</mrow>
</msup>
<mi>&Delta;</mi>
<mi>u</mi>
<msup>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mi>d</mi>
<mi>t</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中,γ1,γ2为任意实数,表示函数f(t)在[t1,t2]上的γ次积分,D为微分符号;
依据Grünwald-Letnikov分数阶微积分定义,对上述目标函数在采样时间TS进行离散化,得到:
J=ZTΛ(γ1,Ts)Z+ΔUTΛ(γ2,Ts)ΔU
其中,
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>&Lambda;</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&gamma;</mi>
<mi>&epsiv;</mi>
</msub>
<mo>,</mo>
<msub>
<mi>T</mi>
<mi>S</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msubsup>
<mi>T</mi>
<mi>S</mi>
<msub>
<mi>&gamma;</mi>
<mi>&epsiv;</mi>
</msub>
</msubsup>
<mi>d</mi>
<mi>i</mi>
<mi>a</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>w</mi>
<mrow>
<mi>P</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>,</mo>
<msub>
<mi>w</mi>
<mrow>
<mi>P</mi>
<mn>2</mn>
</mrow>
</msub>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<msub>
<mi>w</mi>
<mn>1</mn>
</msub>
<mo>,</mo>
<msub>
<mi>w</mi>
<mn>0</mn>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>w</mi>
<mi>q</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>&omega;</mi>
<mi>q</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&gamma;</mi>
<mi>&epsiv;</mi>
</msub>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>&omega;</mi>
<mrow>
<mi>q</mi>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mi>P</mi>
<mo>-</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>&gamma;</mi>
<mi>&epsiv;</mi>
</msub>
<mo>)</mo>
</mrow>
</msubsup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
时,对q<0,ε=1,2;
2.3依据步骤2.2中的目标函数求解得到控制量,形式如下:
ΔU=-(STΛ(γ1,Ts)S+Λ(γ2,Ts))-1SΛ(γ1,Ts)(Gz(k)+ΨΔR)
Δu(k)=[1,0,…,0]ΔU
u(k)=u(k-1)+Δu(k)
2.4在k+l时刻,l=1,2,3,…,依照2.1到2.3中的步骤依次循环求解分数阶模型预测控制器的控制量u(k+l),再将其作用于被控对象。