[发明专利]一种基于正交参数化LTV模型的迭代学习预测控制方法

摘要

本发明公开了一种针对注塑过程保压段基于正交参数化LTV模型的迭代学习预测控制方法。实施步骤如下(1)建立正交参数化LTV模型；(2)模型参数估计；(3)阶次选择；(4)推导ILC‑MPC控制律。本发明充分利用注塑过程重复运行、跟踪轨迹已知的工艺特点，通过引入正交参数化LTV模型，改变了传统建模方法复杂度高、外推性差的特点。在控制过程中，采用迭代学习控制与模型预测控制相结合的控制策略，既可以充分利用历史批次信息，又可以根据预测模型在线进行滚动优化，从而能更迅速更稳定地完成控制需求。

主权项

一种基于正交参数化LTV模型的迭代学习预测控制方法，其特征在于该方法包括以下步骤：步骤(1)、建立正交参数化LTV模型：利用注塑过程重复运行、跟踪轨迹已知的工艺特点,建立正交参数化LTV模型，用以表征注塑过程非线性问题；其中时变系数为时间轴坐标t的非线性函数；步骤(2)、模型参数估计：由于模型精度能在相关控制准则下达到最高，利用Levenberg‑Marquardtt方法求解非线性最小二乘问题确定模型参数；步骤(3)、阶次选择：为了使模型能以最高精度满足控制要求，利用赤池信息准则确定模型阶次；步骤(4)、推导ILC‑MPC控制律：在步骤(1)所获LTV模型基础上求解优化命题，得到最优迭代学习预测控制律，实现注塑机保压段保压压力的控制；步骤(1)具体是针对注塑过程保压段，采用如下LTV模型表示：yk(t)＝G(q,t)uk(t)+vk(t),t＝1,...,N,k＝1,2,...  (1)其中t和k分别代表时间轴坐标和批次轴坐标；N为每个批次时间长度；yk(t),uk(t),vk(t)分别代表第k个批次t时刻的保压压力，阀门开度以及扰动；G(q,t)是从uk(t)到yk(t)的线性时变传递函数；对注塑过程而言，由于扰动存在着很强的批次相关性，因此不可测扰动可用如下沿批次轴的积分白噪声来描述：vk(t)＝vk‑1(t)+wk(t)，其中wk(t)是零均值高斯白噪声；对上述LTV模型即公式(1)，相邻两批次数据做差分，则有yk-(t)=G(q,t)uk-(t)+wk(t)---(2)]]>其中A(q,t)＝1+a1(t)q‑1+...+an(t)q‑n  (4)B(q,t)＝1+b1(t)q‑1+...+bn(t)q‑n  (5)y‾k(t)=yk(t)-yk-1(t)---(6)]]>u‾k(t)=uk(t)-uk-1(t)---(7)]]>上述公式(2)～(7)是LTV‑OE模型结构，其中n为OE模型阶次，ai(t)和bi(t)为时变模型系数；虽然注塑过程具有明显的非线性和时变性，但是其过程存在过程重复运行、跟踪轨迹已知的显著特点，因此ai(t)和bi(t)可以表示为时间轴坐标n的非线性函数，即其中，是一组多项式基函数，m为多项式阶次，和为加权系数；设首项系数为1的Legendre多项式作为正交基函数；得到首项系数为1的Legendre多项式需要经过两步变换；首先，Legendre多项式定义在[‑1,1]上，其一般形式为：P0(x)=1,Pn(x)=12nn!dndxn(x2-1)n---(10)]]>由于Pn(x)中xn的系数为令则可得到首项系数为1的Legendre多项式；其次，由于Legendre多项式定义域为[‑1,1]，而ai(t)和bi(t)中定义在[1,N]上，通过如下变换则可将t∈[1,N]转换为μ∈[‑1,1]；选用首项系数为1的Legendre多项式作基函数，带入公式(8),(9)中得：ai(t)=Σj=1maijP~j-1(t-N+12N-1)---(11)]]>bi(t)=Σj=1mbijP~j-1(t-N+12N-1)---(12);]]>步骤(2)具体是对于LTV‑OE模型,需要估计的模型参数为：θ=[a11,a12,...,a1m,...,an1,an2,...,anm,b11,b12,...,b1m,...,bn1,bn2,...,bnm]T---(13)]]>则一步向前最优预报值为：y‾^k(t|θ)=B(q,t)A(q,t)u‾k(t)---(14)]]>因此，通过最小化预报误差准则函数VKN(θ)，来确定参数矢量θ：VKN(θ)=1KΣk=1K[1NΣt=1Nϵk(t|θ)2]---(15)]]>其中，由于εk(t|θ)与A(q,t)呈非线性关系，则优化命题(15)不存在解析解，需采用数值方法进行求解；对于上述非线性最小二乘问题，采用L‑M法进行数值求解，并得到如下迭代计算过程：θ^r+1=θ^r+[Σk=1KΣt=1Nψk(t,θ^r)ψkT(t,θ^r)+μI]-1[Σk=1KΣt=1Nψk(t,θ^r)ϵk(t|θ^r)]---(16)]]>其中为第r次迭代的估计值，μ为步长因子，I为单位阵，ψk(t,θ)=-[∂∂θϵk(t|θ)]T=[∂∂θy‾^k(t|θ)]T---(17)]]>在用上述方法进行模型参数的数值优化时，需知道的梯度ψk(t,θ)；对于公式(14)，两边分别对和求导得：A(q,t)·∂y‾^k(t|θ)∂aij=-P~j-1(t-N+12N-1)·y‾^k(t-i|θ)---(18)]]>A(q,t)·∂y‾^k(t|θ)∂bij=P~j-1(t-N+12N-1)·u‾k(t-i)---(19)]]>将A(q,t)代入上式，可得梯度ψk(t,θ)的计算公式；步骤(3)具体是采用赤池信息准则确定LTV‑OE模型的模型阶次n以及正交多项式阶次m；所述的赤池信息准则为：AIC=logVKN(θ)+2dNk---(20)]]>其中，Nk为辨识数据个数，Nk＝K·N；d为模型参数个数，d＝2n·m；步骤(4)具体是对于LTV模型相邻批次做差分有：yk(t)-yk-1(t)=B(q,t)A(q,t)[uk(t)-uk-1(t)]+vk(t)-vk-1(t)---(21)]]>即yk(t)-yk-1(t)=B(q,t)A(q,t)rk(t)+wk(t)---(22)]]>写成状态空间形式为：xk(t+1)=A(t+1)xk(t)+B(t+1)rk(t)yk(t)-yk-1(t)=Cxk(t)+wk(t)---(23)]]>其中，预测模型：y^k(t+j|t)-yk-1(t+j)=Cxk(t+j|t)+w^k(t+j|t)---(25)]]>令则y^k(t+j|t)-yk-1(t+j)=Cxk(t+j|t)+wk(t)---(26)]]>yk(t)‑yk‑1(t)＝Cxk(t)+wk(t)因此有：y^k(|t+Pt+1,t)=yk-1(|t+Pt+1)+Ipx1·C·xk(|t+Pt+1,t)+Ipx1·wk(t)---(27)]]>P为预测时域；由于xk(t)为已知状态，则根据递推公式有：xk(t+1|t)=A(t+1)xk(t)+B(t+1)r^(t)---(28)]]>xk(t+2|t)=A(t+2)xk(t+1|t)+B(t+2)rk(t+1|t)=A(t+2)A(t+1)xk(t)+A(t+2)B(t+1)r^k(t)+B(t+2)r^k(t+1|t)...xk(t+M|t)=A(t+M)...A(t+1)xk(t)+A(t+M)...A(t+2)B(t+2)r^k(t)+...+B(t+M)r^k(t+M-1|t)xk(t+M+1|t)=A(t+M+1)...A(t+1)xk(t)+A(t+M+1)...A(t+2)B(t+1)r^k(t)+...+A(t+M+1)B(t+M)r^k(t+M-1|t)...xk(t+P|t)=A(t+P)...A(t+1)xk(t)+A(t+P)...A(t+2)B(t+1)r^k(t)+...+A(t+P)...A(t+M+1)B(t+M)r^k(t+M-1|t)---(29)]]>由于控制时域为M，则有：综上推导，可得，xk(t+1|t)...xk(t+M|t)xk(t+M+1|t)...xk(t+P|t)=A(t+1)...Πi=1MA(t+i)Πi=1M+1A(t+i)...Πi=1PA(t+i)xk(t)+B(t+1)0...0............Πi=2MA(t+i)B(t+1)......B(t+M)Πi=2M+1A(t+i)B(t+1)......A(t+M+1)B(t+M)............Πi=2PA(t+i)B(t+1)......Πi=M+1PA(t+i)B(t+M)r^k(t)............r^k(t+M-1|t)]]>令h(|t+Pt+1)A(t+1)...Πi=1MA(t+i)Πi=1M+1A(t+i)...Πi=1PA(t+i),g(|t+Pt+1)=B(t+1)0...0............Πi=2MA(t+i)B(t+1)......B(t+M)Πi=2M+1A(t+i)B(t+1)......A(t+M+1)B(t+M)............Πi=2PA(t+i)B(t+1)......Πi=M+1PA(t+i)B(t+M)]]>则可表示为：xk(|t+Pt+1,t)=h(|t+Pt+1)xk(t)+g(|t+Pt+1)r^k(|t+M-1t,t)---(30)]]>带入(27)有：y^k(|t+Pt+1,t)=yk-1(|t+Pt+1)+Ipx1·C·[h(|t+Pt+1)·xk(t)+g(|t+Pt+1)·r^k(|t+M-1t,t)]+IPx1·wk(t)]]>整理得，y^k(|t+Pt+1,t)=G(|t+Pt+1)r^k(|t+M-1t,t)+fk(t+Pt+1)---(31)]]>其中，G(|t+Pt+1)=CB(t+1)0...0............CΠi=2MA(t+i)B(t+1)......CB(t+M)CΠi=2M+1A(t+i)B(t+1)......CA(t+M+1)B(t+M)............CΠi=2PA(t+i)B(t+1)......CΠi=M+1PA(t+i)B(t+M)]]>fk(t+Pt+1)=yk-1(|t+Pt+1)+H(|t+Pt+1)·xk(t)+IPx1·wk(t)---(32)]]>选目标函数为：J(t,k,M,P)=||yr(|t+Pt+1)-y^k(|t+Pt+1,t)||Q2+||rk(|t+M-1t)||R2---(33)]]>无约束问题解析解为：rk*=(GTQG+R)-1GTQ·[yr(|t+Pt+1)-f(|t+Pt+1)]---(34)]]>令dT＝[1 0 … 0](GTQG+R)‑1GTQ,则rk=dT[yr(|t+Pt+1)-f(|t+Pt+1)]---(35)]]>fk(t+Pt+1)=yk-1(|t+Pt+1)+H(|t+Pt+1)·xk(t)+IPx1·wk(t)---(36)]]>其中，H(|t+Pt+1)=CA(t+1)...CΠi=1MA(t+i)CΠi=1M+1A(t+i)...CΠi=1PA(t+i),G(|t+Pt+1)=CB(t+1)0...0............CΠi=2MA(t+i)B(t+1)......CB(t+M)CΠi=2M+1A(t+i)B(t+1)......CA(t+M+1)B(t+M)............CΠi=2PA(t+i)B(t+1)......CΠi=M+1PA(t+i)B(t+M)]]>求解优化命题公式(33)，得到最终解析解公式(35)，从而得到控制律为：uk(t)＝rk(t),k＝1；uk(t)＝rk(t)+uk‑1(t),k＞1  (37)。