CN103995466A - 一种基于软约束的区间预测控制建模及优化方法 - Google Patents

Abstract

一种基于软约束的区间预测控制建模及优化方法。该控制方法步骤如下：(1)基于过程预测模型，建立包含约束项、控制项和经济项的二次性能指标；(2)通过求解松弛变量，判断综合优化方法是否可行(3)给出了控制模型输出约束不可行时软约束松弛变量的求解方法，实现了区间预测控制模型输出约束不可行时可行域范围的调整；(4)采用一种边界可行序列二次规划方法，用以解决初始点选择不佳导致方法增加计算量或难以找到最优解以及计算中舍入误差的影响会破坏Hessian矩阵正定性等情形的问题，并求出最优控制输入。本发明能够建立复杂的多变量系统控制模型，基于软约束调整准确快速求解出控制律，有利于实现多变量系统的良好控制。

Description

一种基于软约束的区间预测控制建模及优化方法

技术领域

本发明属于过程控制领域，尤其涉及一种基于软约束的区间预测控制建模及优化方法。

背景技术

近年来，随着工业过程系统的规模日趋庞大和能源供应的日趋紧张，以及能源价格一再拉高，生产者们不再单单对生产过程的某一参量或某一性能提出控制要求，而是根据生产要求以及人为意愿提出经济性、快速性等综合性能指标控制。然而由于环境因素和系统自身的复杂性，实际系统存在多种约束，因而普通的预测控制方法是很难满足控制要求。

在实际工业控制过程中，如果把系统区域约束条件考虑到控制目标中会导致优化控制器的可行集缩小，或者多层优化指标的引入又会使控制器求解计算量很大，所以传统的优化方法不能同时考虑优化经济和处理约束目标。考虑到控制模型算法的可行性，早些时候的线性化处理方法，凸多面体几何角度方法都是针对线性模型，对在非线性系统上的应用是否可行尚未得到验证。另外，对于求解算法，各种改进形式的遗传算法、粒子群算法等多种非线性优化求解算法。这些算法在求解控制模型时需要多次迭代计算，因此计算量很大，程序运行起来繁琐复杂，从而导致实际系统快速运行时无法实现实时控制。

综上所述，现有方法不是有局限性就是求解计算量大，很难应用到实践中去。本发明对于复杂的多变量有约束系统，建立了一个能够处理约束目标、协调操作控制和改善经济指标控制模型，并且设计了一个既能够自动放松输出约束可行集又能快速、准确求解控制模型的方法。

发明内容

为了克服上述现有技术的不足，本发明提供了一种可以解决经济性和控制性协调优化的控制模型，当约束集不可行时能快速准确求解控制输入方法。

本发明解决其技术问题所采用的技术方案是：

一种基于软约束的区间预测控制建模及优化方法，其内容包括如下步骤：

步骤一：利用被控变量(controlled variable，缩写为CV)预测输出与软、硬约束构造区间约束项指标，对输出进行区间控制；

步骤二：记CV预测输出超出软约束区间的上线或下线的偏差值为松弛变量A1或A2，软约束上线与硬约束上线和软约束下线与硬约束下线差值分别记为B1和B2；根据A1与B1或A2与B2的比值确定区间约束项加权矩阵元素的大小：当比值大于1时，对应的权值矩阵元素都为1；当比值介于0和1之间时，对应的权值矩阵元素为比值的平方；当比值为0时，对应的权值矩阵元素都为0；

步骤三：利用CV经济约束构造CV区间控制项指标，利用操作变量(manipulated variable，缩写为MV)经济约束构造MV区间控制项指标，CV区间控制项指标与MV区间控制项指标之和统称为区间控制项指标；

步骤四：利用输出量、控制量和控制量增量的二次性能指标之和构造控制模型经济项指标；将所述的输出量和控制量的权值矩阵元素与产品产量有关的权值设为负值，将与成本有关的权值设为正值，而将与产量和成本均无关的二次指标对应的权值设为零，则控制增量的加权矩阵为单位矩阵；

步骤五：用区间约束项指标、区间控制项指标和经济项指标的二次性能指标三者之和作为综合优化方法的控制模型，各项权值矩阵的调节遵循变量优先权原则，优先级顺序为：约束项、经济项、控制项；只有被控变量处于软约束范围内时才能对经济项和控制项的权重进行协调，此时区间约束项权值全部为0，且经济项的权值的最小值大于对应项控制项的权值的最大值；

步骤六：通过求解松弛变量A1或A2是否为0，判断求解方法是否可行。当求解方法不可行时，输出软约束放松大小为松弛变量A1或A2，使控制模型重新变得可行；在可行性分析后，当松弛变量的解为0时，约束控制模型可行，当松弛变量的解不为0时，约束控制模型不可行；

步骤七：当求解方法可行时，若初始点选择不佳，用蒙特卡罗方法找到接近最优点的可行域点；

步骤八：在Hessian矩阵正定性的确保下，采用积极约束估计集方法，减少二次规划方法中子问题的约束条件，并结合最速下降与拟牛顿方法求取最优化过程的可行下降方向，通过循环迭代求取最优控制输入。

在步骤七中所述的蒙特卡罗方法，记为方法一，其内容包括如下步骤：

第一步预置N为充分大的正数，确定选点个数M；

第二步用随机函数及条件限制产生可行点x；

第三步计算控制模型:F＝f(x)；

第四步比较函数值:若F≥N，转第六步；否则，转第五步；

第五步记录当前最优点的信息:N＝F,x^k＝x；

第六步若已选定M个可行点，输出x^k和N；否则，转第二步，寻找下一个可行点。

在上述步骤八中所述的积极约束估计集方法，记为方法二，其内容包括如下步骤：

步骤1对于x^k∈X，令i＝0,ε_k,i(x^k)＝ε₀∈(0,1)。

步骤2令

$\begin{array}{c}{L}_{k,i}\left(x^k\right)=\{j\∈I|-{\mathrm{\ϵ}}_{k,i}\left(x^k\right)\≤g_j\left(x^k\right)\≤0,{\mathrm{\λ}}_j\≥0\},\\ A_{k,i}\left(x^k\right)={\▿g}_{L_{k,i}\left(x^k\right)}\left(x^k\right)\end{array}$

步骤3如果det(A_k,i(x^k)^TA_k,i(x^k))<ε_k,i(x^k)，其中，

L_k＝φ，A_k＝φ，i_k＝φ，i＝i+1，ε_i(x^k)＝ε_i(x^k)/2，

返回步骤2；否则，令：L_k＝L_k,i(x^k)，A_k＝A_k,i(x^k)，i_k＝i，则循环停止。

由于此方法中获得的A_k为列满秩矩阵,寻找其行的最大线性无关组,并对其行进行重新排列,记:

$A_k=\left(\begin{array}{c}{A}_k^1\\ A_k^2\end{array}\right),\&dtri;f\left(x^k\right)=\left(\begin{array}{c}\&dtri;f^1\left(x^k\right)\\ {\&dtri;f}^2\left(x^k\right)\end{array}\right),$

其中是可逆矩阵,并且是A_k的|L_k|个线性独立的行合成的矩阵,那么为A_k剩余n-|L_k|行合成的矩阵。类似A_k,f(x^k)的组成同理。

在所述的步骤八中采用重置H_k为单位矩阵保证Hessian矩阵正定性，此时的迭代点作为初始点的方法重新进行迭代.具体更新H_k的改进BFGS公式如下，记为方法三：

$H_{k+1}=\left\{\begin{array}{cc}I,& \left|\right|x^{k+1}-x^k\left|\right|<\mathrm{\&epsiv;}\\ H_k+\frac{\widetilde{y_k}\widetilde{{y_k}^T}}{\widetilde{{y_k}^T}{S}_k}-\frac{H_k{s}_k{s}_k^T{H}_k}{s_k^T{H}_k{s}_k},& \left|\right|x^{k+1}-x^k\left|\right|\&GreaterEqual;\mathrm{\&epsiv;}\end{array}\right.---\left(1\right)$

其中,ε为一个给定的极小值,

$s_k=x_{k+1}-x_k,y_k={\&dtri;}_{x}L(x_{k+1},{\mathrm{\&lambda;}}_{k+1})-{\&dtri;}_{x}L(x_k,{\mathrm{\&lambda;}}_{k+1})$

$\mathrm{\&theta;}=\left\{\begin{array}{cc}1& (s_k^T{y}^k\&GreaterEqual;0.2s_k^T{H}_s{s}_k)\\ \frac{0.8s_k^T{H}_k{s}_k}{s_k^T{H}_k{s}_k-s_k^T{y}_k}& (s_k^T{y}^k<0.2s_k^T{H}_k{s}_k)\end{array}\right.---\left(2\right)$

$\widetilde{y_k}=\mathrm{\&theta;}{y}_k+(1+\mathrm{\&theta;})H_k{s}_k$

本发明的有益效果主要表现在：

1.控制模型由约束项、控制项和经济项构成，不仅能够处理约束目标，而且可以改善经济指标，能实现多变量控制系统的经济性能优化以及安全顺利运行；

2.根据被控变量超出软硬约束的程度，控制模型约束项权值可以自动调整，使被控变量可以很快调整到约束范围内，保证了产品质量；

3.在软约束下当控制器求解不可行时，软约束范围可以根据被控变量需要自动调整，保证控制器的正常工作；

4.设计了边界可行序列二次规划求解方法。该方法计算量小、计算速度快、计算精度高，能够准确求解带边界约束的非线性规划最优化问题，从而实现了区间预测控制模型的快速、精确求解，对实现多变量系统的良好控制具有重要的实用价值。

附图说明

图1为CV区间约束项的构造示意图；

图2为CV区间控制项的构造示意图；

图3为CV软约束调整后约束区间构造示意图。

具体实施例

下面结合附图和具体实施例对本发明做进一步的详细描述。

以壳牌石油重油分馏塔的三入三出系统为例。过程的MV：u₁代表分馏器顶部产品的抽出率；u₂代表分馏器侧线产品的抽出率；u₃代表分馏器底部的回流热负荷。过程的CV：y₁代表分馏器顶部产品的提取成分；y₂代表分馏器侧线产品的提取成分；y₃代表分馏器底部的回流温度。取预测步数为20，控制步数为10。

1、构造对重油分馏塔的石油提取速率的进行约束的约束项指标J₁

CV区间约束项构造如图1所示，一般情况下，点7在线2与线3之间便不再受惩罚。当分馏塔中部和塔顶的提取成分过高时，CV预测值超出线1或线4之后则被高度惩罚；当点6介于线1与线2之间，或当点8介于线3与线4之间时，CV的惩罚程度有所降低。粗实线可以表示松弛变量的值，用字母和表示，公式如下：

$S_C^{+}=\widehat{Y_c}\left(k\right)-Y_{\mathrm{HH}}---\left(1\right)$

$S_C^{-}=Y_{\mathrm{LL}}-\widehat{Y_c}\left(k\right)---\left(2\right)$

分别表示对应和的权值对角矩阵，且

$P_C^{+}=\mathrm{diag}\{p_{C1}^{+},p_{C2}^{+},\&CenterDot;\&CenterDot;\&CenterDot;,p_{C20}^{+}\},p_{\mathrm{Cj}}^{+}=\mathrm{diag}\{p_{\mathrm{Cj}1}^{+},p_{\mathrm{Cj}2}^{+},p_{\mathrm{Cj}3}^{+}\},(j=\mathrm{1,2}\&CenterDot;\&CenterDot;\&CenterDot;,20)$

$P_C^{-}=\mathrm{diag}\{p_{C1}^{-},p_{C2,}^{-}\&CenterDot;\&CenterDot;\&CenterDot;,p_{C20}^{-}\},p_{\mathrm{Cj}}^{-}=\mathrm{diag}\{p_{\mathrm{Cj}1}^{-},p_{\mathrm{Cj}2}^{-},p_{\mathrm{Cj}3}^{-}\},(j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,20)$

则约束项指标为下式：

$J_1=\left|\right|S_C^{-}|{|}_{P_C^{-}}^2+|\left|S_C^{+}\right|{|}_{P_C^{+}}^2---\left(3\right)$

2、根据偏差值自动调节约束项权值矩阵，保证分馏器抽出率平稳

本发明定义松弛变量为预测输出超出(高于或者低于)约束区间的数值，全部为非负实数。表示校正后的被控变量预测值。

对于对角阵和中的元素值，根据偏差值来确定它们的大小。首先定义如下四个偏差值变量：

${\mathrm{\&Delta;}}_{\mathrm{HHij}}=\widehat{Y_{\mathrm{ci}}}(k+j|k)-Y_{\mathrm{HHi}}\left(j\right)$

θ_HHij＝Y_HHHi(_j)-Y_HHi(_j)   (4)

${\mathrm{\&Delta;}}_{\mathrm{LLij}}=Y_{\mathrm{LLi}}\left(j\right)-\widehat{Y_{\mathrm{ci}}}(k+j|k)$

θ_LLij＝Y_LLi(j)-Y_LLLi(j)

式中，Δ_HHij表示第i个输出的第j步预测值超出软约束上限的偏差值，θ_HHij表示第i个输出的第j步预测值对应的硬约束上线与软约束上限的偏差值，Δ_LLi_j表示第i个输出的第j步预测值超出软约束下限的偏差值，θ_LLij表示第i个输出的第j步预测值对应的硬约束下线与软约束下限的偏差值。其中各项元素的值具体设定如下：

另x＝Δ_HHij/θ_HHij,y＝Δ_LLij/θ_LLij

$f\left(x\right)=\left\{\begin{array}{cc}{x}^2,& 0<x<1\\ 1,& x\&GreaterEqual;1\end{array}\right.---\left(5\right)$

$f\left(y\right)=\left\{\begin{array}{cc}{y}^2,& 0<y<1\\ 1,& y\&GreaterEqual;1\end{array}\right.---\left(6\right)$

如图1所示：

(1)若点5超出硬约束上线1，则Δ_HHij/θ_HHij≥1，那么

(2)若点6介于软约束上线2与硬约束上线1之间，则0<Δ_HHij/θ_HHij<1，那么 $P_{\mathrm{Cji}}^{-}=0;$

(3)若点24介于软约束上下线2与3之间，则和均为0。即在此区间中，此时输出在此区间内的变化就可忽略。

(4)同理，若预测输出处于线3与线4之间，

(5)同理，若预测输出超出线4，

这样就能够保证每次求解式(3)描述的区间控制模型时，其各项权值均能根据偏差大小自适应调节。

3、构造区间控制项指标构造J₂，协调操作石油提取速率和塔底回流温度，以保证较好的经济效益

多变量系统控制模型的区间控制项，可分为CV区间控制项与MV区间控制项。

对于式(6)描述的控制模型区间约束项，CV在线1与线4之间(包括边界)后便不再受惩罚。为使分馏塔提取石油等产品产量更大，这就需要被控变量被控制在线8与线9的小区间内，如图2所示。CV在线8与线9之间不受惩罚，超出此范围接收惩罚。如图2所示，粗实线为松弛变量，字母表示为P_Cbelow和P_Cabove。对于每一个预测输出，其双限小区间控制的等式形式如下：

$S_{\mathrm{Cbelow}}^{-}=Y_L-\widehat{Y_c}\left(k\right)---\left(7\right)$

$S_{\mathrm{Cabove}}^{+}=\widehat{Y_c}\left(k\right)-Y_H---\left(8\right)$

式中，Y_H和Y_L分别表示为线8与线9，是双限小区间的上下界。因此，应被加到控制模型中，其中P_Cbelow,P_Cabove∈R^(3□20)×(3^□20)分别表示对应和的权值对角矩阵，且

P_Cbelow＝diag{p_cbelow1,p_cbelow2,…,p_cbelow20},p_cbelowj＝diag{p_cbelowj1,p_cbelowj2,p_cbelowj3},

P_cabove＝diag{p_cabove1,p_cabove2,…,p_cabove20},p_cabovej＝diag{p_cabovej1,p_cabovej2,p_cabovej3},(j＝1,2,…,20).

惩罚松弛变量二次方方法在建立CV区间控制项中的应用对于实现MV的区间控制(处理MV区间控制项)也同样适用。因此，在MV约束范围内定义软约束，同理可在控制模型中加入MV区间控制项

D_Mabove＝diag{d_mabove1,d_mabove2,…,d_mabove10},d_mabovei＝diag{d_mabovei1,d_mabovei2,d_mabovei3},

D_Mbelow＝diag{d_mbelow1,d_mbelow2,…,d_mbelow10},d_mbelowi＝diag{d_mbelowi1,d_mbelowi2,d_mbelowi3},(i＝1,2,L,10)

可得控制模型的区间控制项公式为：

$J_2=\left|\right|S_{\mathrm{Cbelow}}^{-}|{|}_{P_{\mathrm{Cbelow}}}^2+|\left|S_{\mathrm{Cabove}}^{+}\right|{|}_{P_{\mathrm{Cabove}}}^2+\left|\right|S_{\mathrm{Mbelow}}^{-}|{|}_{D_{\mathrm{Mbelow}}}^2+|\left|S_{\mathrm{Mabove}}^{+}\right|{|}_{D_{\mathrm{Mabove}}}^2---\left(9\right)$

4、构造产品成本与产量的经济项，达到经济效益最大化

为满足对生产能耗与经济收益的优化，本节构造了如下用于协调石油提取成分与抽出率的多变量系统控制模型的经济项：

$J_3=\left|\right|\widehat{Y_c}\left(k\right)|{|}_{P_C}^2+|\left|U\left(k\right)\right|{|}_D^2+\left|\right|\mathrm{\&Delta;U}\left(k\right)|{|}_{\mathrm{DD}}^2---\left(10\right)$

式中，P_C∈R^(3□20)×(3^□20),D,DD∈R^{(3□10)×(3□10)}分别表示二次性能指标||U(k)||²、||ΔU(k)||²对应的权值矩阵，它们均为对角矩阵，用来对这些二次指标施以惩罚，且

P_c＝diag{p_c1,p_c2,…,p_c20},p_cj＝diag{p_cj1,p_cj2,…,p_cj3},(j＝1,2,…,20),

p_c(j+1)＝10^-2p_cj,(j＝1,2,…19)；

D＝diag{d₁,d₂,…,d₁₀},d_i＝diag{d_i1,d_i2,d_i3},(i＝1,2,…,10),

d_i+1＝10^-2d_i,(i＝1,2,…,9)；

DD＝diag{dd₁,dd₂,dd₁₀},dd_i＝diag{dd_i1,dd_i2,dd_i3},(i＝1,2,…,10).

权值对角矩阵P_C和D中，与分馏塔产品产量有关的权值设为负值，以保证其负收益最小，即正收益最大；与分馏塔控制成本有关的权值设为正值，以保证其消耗成本最低，且这些权值的绝对值均大于控制项中各权值，而与产量和成本均无关的二次指标对应的权值设为零。通过调节P_C和D可在约束区间内协调优化CV与MV值，从而使经济收益实现一定程度的提高。

由于模型预测输出与输入之间存在着一定的增益关系，这就意味着所有的输出变量增量均可由某些或全部的输入变量增量表示，因此式(10)描述的控制模型经济项中用来限制MV在系统运行过程中的动态行为。

5、预测控制模型构造

该控制模型考虑3个性能指标，组成形式如下：

J＝J₁+J₂+J₃    (11)

式中，J为控制模型的总体性能指标，J₁表示区间约束项指标，表示对重油分馏塔的石油提取速率的约束；J₂表示区间控制项指标，为保证较好的经济状态，协调操作石油提取速率和塔底回流温度；J₃表示区间经济项指标，确保石油提取率和回流温度达到最大化经济值。

综上所述，将式(3)、(9)和(10)代入式(12)，可得下列多变量系统的区间预测控制模型：

$\begin{array}{c}J\left(k\right)=\left|\right|\widehat{Y_c}\left(k\right)-Y_{\mathrm{HH}}|{|}_{P_C^{+}}^2+||Y_{\mathrm{LL}}-\widehat{Y_c}\left(k\right)|{|}_{P_C^{-}}^2+\left|\right|\widehat{Y_c}\left(k\right)-Y_H|{|}_{P_{\mathrm{Cabove}}}^2+\\ \left|\right|Y_L-\widehat{Y_c}\left(k\right)|{|}_{P_{\mathrm{Cbelow}}}^2+||U\left(k\right)-U_H|{|}_{D_{\mathrm{Mabove}}}^2+\left|\right|U_L-U\left(k\right)|{|}_{D_{\mathrm{Mbelow}}}^2+\\ \left|\right|\widehat{Y_c}\left(k\right)|{|}_{P_C}^2+|\left|U\left(k\right)\right|{|}_D^2+\left|\right|\mathrm{\&Delta;U}\left(k\right)|{|}_{\mathrm{DD}}^2\end{array}a---\left(12\right)$

s.t.ΔU_LL≤ΔU(k)≤ΔU_HH

6、控制方法可行性研究

通过图3所示，为保证式(12)描述的约束控制模型可行，对输出约束引入松弛变量ε₁和ε₂，所以原来的软约束上线2变为图3中的线13，软约束下线3变为图3中线14。通过QP方法最小化其中若所得最优解ε＝0，则式(12)描述的控制模型约束可行；若所得最优解ε≠0，则说明式(12)描述的控制模型约束不可行，此时需对软约束放松ε，使控制模型重新变的可行。

引入松弛变量对控制模型进行改进，新的软约束上下线分别为和

$Y_{\mathrm{HH}}^0=Y_{\mathrm{HH}}+{\mathrm{\&epsiv;}}_1---\left(13\right)$

$Y_{\mathrm{LL}}^0=Y_{\mathrm{LL}}-{\mathrm{\&epsiv;}}_2---\left(14\right)$

将式(12)中约束项中的Y_HH替换为Y_LL替换为对偏差进行惩罚，可用和两项来实现。可得到如下能够进行输出软约束调整的约束控制模型：

$\begin{array}{c}J\left(k\right)=\left|\right|\widehat{Y_c}\left(k\right)-Y_{\mathrm{HH}}|{|}_{P_C^{+}}^2+||Y_{\mathrm{LL}}-\widehat{Y_c}\left(k\right)|{|}_{P_C^{-}}^2+\left|\right|\widehat{Y_c}\left(k\right)-Y_H|{|}_{P_{\mathrm{Cabove}}}^2+\\ \left|\right|Y_L-\widehat{Y_c}\left(k\right)|{|}_{P_{\mathrm{Cbelow}}}^2+||U\left(k\right)-U_H|{|}_{D_{\mathrm{Mabove}}}^2+\left|\right|U_L-U\left(k\right)|{|}_{D_{\mathrm{Mbelow}}}^2+\\ \left|\right|\widehat{Y_c}\left(k\right)|{|}_{P_C}^2+|\left|U\left(k\right)\right|{|}_D^2+\left|\right|\mathrm{\&Delta;U}\left(k\right)|{|}_{\mathrm{DD}}^2\end{array}a---\left(15\right)$

s.t.ΔU_L≤ΔU(k+i)≤ΔU_H

7、选择最佳初始点

由蒙特卡罗方法(方法一)寻找一个比较好的初始点U⁰∈R³；

第一步设N为1.0000e+010，确定选点个数为10；

第二步用随机函数和模型约束条件产生可行点U；

第三步计算约束控制模型F＝J(k)；

第四步比较函数值:若F≥N，转第六步；否则，转第五步；

第五步记录当前最优点的信息:N＝F,U^k＝U；

第六步若10个可行点已经选出，则输出U^k和N，否则，转第二步，寻找下一个可行点。

8、结合最速下降与拟牛顿方法寻优，通过循环迭代求取最优控制输入

针对式(15)描述的区间预测控制模型的特点，提出如下边界可行序列二次规划方法(BFSQP)，其内容包括如下步骤：

步骤1)由方法二获得近似积极估计集；所述步骤八中所述的积极约束估计集方法；

步骤2)k＝1时，令步长t_k＝1，若t_k>γ,则的求取采用最速下降法,即令否则通过求解小规模QP子问题获得KKT点对如果转步骤5,否则,停止；

步骤3)选择合适的M₁(U^k),M₂(U^k)和M₃(U^k)表达式,并通过 $p^k=-M_1\left(U^k\right){\left({\left(A_k^1\right)}^{-1}\right)}^{T}e,d_1^k=\left(\begin{array}{c}{p}^k\\ 0\end{array}\right),d^k=M_2\left(U^k\right)d_0^k+M_3\left(U^k\right)d_1^k$ 获得可行下降方向d^k；

步骤4)由 ${\left(A_k^1\right)}^T\mathrm{\&rho;}=-\left|\right|d_0^k|{|}^{\mathrm{\&tau;}}e-\widehat{g_{L_k}}(x^k+d^k)$ 和 ${\left(A_k\right)}^T\widehat{d^k}={\left(A_k^1\right)}^T{\mathrm{\&rho;}}^k+{\left(A_k^2\right)}^{T}0={\left(A_k^1\right)}^T{\mathrm{\&rho;}}^k$ 获得一个高阶修正方向

步骤5)计算最大值使其满足以下两条件:

$\begin{array}{c}f(U^k+t_k{d}^k+t_k^2\widehat{d^k})\&le;f\left(U^k\right)+\mathrm{\&alpha;}{t}_k\&dtri;f\left(U^k\right)d^k,\\ g_j(U^k+t_k{d}^k+t_k^2\widehat{d^k})\&le;0,j\&Element;I\end{array}$

步骤6)迭代:

λ^k+1＝λ^k+t_k△λ^k；

${\mathrm{\&mu;}}^{k+1}\left(i\right)=\min \left\{\max \right\{{\mathrm{\&lambda;}}^{k+1}\left(i\right),{\mathrm{\&lambda;}}^{k0}\left(i\right),\left|\right|d^k\left|\right|\},\stackrel{\&OverBar;}{\mathrm{\&mu;}}\},i\&Element;L_k,$

μ^k+1(i)＝μ^k(i)/2,i∈I/L_k

若积极估计集L_k≠φ，则否则U^k+1＝U^k+t_kd^k；

步骤7)根据方法三更新对称正定矩阵H_k获得H_k+1,k＝k+1,转入步骤2)，直至迭代停止。

Claims (4)

1.一种基于软约束的区间预测控制建模及优化方法，其特征在于：该方法内容包括如下步骤：

步骤一：利用被控变量预测输出与软、硬约束构造区间约束项指标，对输出进行区间控制；

步骤二：记被控变量预测输出超出软约束区间的上线或下线的偏差值为松弛变量A1或A2，软约束上线与硬约束上线和软约束下线与硬约束下线差值分别记为B1和B2；根据A1与B1或A2与B2的比值确定区间约束项加权矩阵元素的大小：当比值大于1时，对应的权值矩阵元素都为1；当比值介于0和1之间时，对应的权值矩阵元素为比值的平方；当比值为0时，对应的权值矩阵元素都为0；

步骤三：利用被控变量经济约束构造被控变量区间控制项指标，利用操作变量经济约束构造操作变量区间控制项指标，被控变量区间控制项指标与操作变量区间控制项指标之和统称为区间控制项指标；

步骤四：利用输出量、控制量和控制量增量的二次性能指标之和构造控制模型经济项指标；将所述的输出量和控制量的权值矩阵元素与产品产量有关的权值设为负值，将与成本有关的权值设为正值，而将与产量和成本均无关的二次指标对应的权值设为零，则控制增量的加权矩阵为单位矩阵；

步骤五：用区间约束项指标、区间控制项指标和经济项指标的二次性能指标三者之和作为综合优化方法的控制模型，各项权值矩阵的调节遵循变量优先权原则，优先级顺序为：约束项、经济项、控制项；只有被控变量处于软约束范围内时才能对经济项和控制项的权重进行协调，此时区间约束项权值全部为0，且经济项的权值的最小值大于对应项控制项的权值的最大值；

步骤六：通过求解松弛变量A1或A2是否为0，判断求解方法是否可行；当求解方法不可行时，输出软约束放松大小为松弛变量A1或A2，使控制模型重新变得可行；在可行性分析后，当松弛变量的解为0时，约束控制模型可行，当松弛变量的解不为0时，约束控制模型不可行；

步骤七：当求解方法可行时，若初始点选择不佳，用蒙特卡罗方法找到接近最优点的可行域点；

步骤八：在Hessian矩阵正定性的确保下，采用积极约束估计集方法，减少二次规划方法中子问题的约束条件，并结合最速下降与拟牛顿方法求取最优化过程的可行下降方向，通过循环迭代求取最优控制输入。

2.根据权利要求1所述的一种基于软约束的区间预测控制建模及优化方法，其特征在于：在所述步骤七中所述的蒙特卡罗方法，其内容包括如下步骤：

第一步预置N为充分大的正数，确定选点个数M；

第二步用随机函数及条件限制产生可行点x；

第三步计算控制模型:F＝f(x)；

第四步比较函数值:若F≥N，转第六步；否则，转第五步；

第五步记录当前最优点的信息:N＝F,x^k＝x；

第六步若已选定M个可行点，输出x^k和N；否则，转第二步，寻找下一个可行点。

3.根据权利要求1所述的一种基于软约束的区间预测控制建模及优化方法，其特征在于：在所述步骤八中所述的积极约束估计集方法，其内容包括如下步骤：

步骤1对于x^k∈X，令i＝0,ε_k,i(x^k)＝ε₀∈(0,1)，

步骤2令

$\begin{array}{c}{L}_{k,i}\left(x^k\right)=\{j\&Element;I|-{\mathrm{\&epsiv;}}_{k,i}\left(x^k\right)\&le;g_j\left(x^k\right)\&le;0,{\mathrm{\&lambda;}}_j\&GreaterEqual;0\},\\ A_{k,i}\left(x^k\right)={\&dtri;g}_{L_{k,i}\left(x^k\right)}\left(x^k\right)\end{array}$

步骤3如果det(A_k,i(x^k)^TA_k,i(x^k))<ε_k,i(x^k)，其中，

L_k＝φ，A_k＝φ，i_k＝φ，i＝i+1，ε_i(x^k)＝ε_i(x^k)/2，

返回步骤2；否则，令：L_k＝L_k,i(x^k)，A_k＝A_k,i(x^k)，i_k＝i，则循环停止；

由于此方法中获得的A_k为列满秩矩阵,寻找其行的最大线性无关组,并对其行进行重新排列,记:

$A_k=\left(\begin{array}{c}{A}_k^1\\ A_k^2\end{array}\right),\&dtri;f\left(x^k\right)=\left(\begin{array}{c}\&dtri;f^1\left(x^k\right)\\ {\&dtri;f}^2\left(x^k\right)\end{array}\right),$

其中是可逆矩阵,并且是A_k的|L_k|个线性独立的行合成的矩阵,那么为A_k剩余n-|L_k|行合成的矩阵；类似A_k,f(x^k)的组成同理。

4.根据权利要求1所述的一种基于软约束的区间预测控制建模及优化方法，其特征在于：在所述步骤八中所述的二次规划方法，其内容包括如下步骤：

步骤1)由所述步骤八中所述的积极约束估计集方法获得近似积极估计集；

步骤2)k＝1时，令步长t_k＝1，若t_k>γ,则的求取采用最速下降法,即令否则通过求解小规模QP子问题获得KKT点对如果转步骤5,否则,停止；

步骤3)选择合适的M₁(U^k),M₂(U^k)和M₃(U^k)表达式,并通过 $p^k=-M_1\left(U^k\right){\left({\left(A_k^1\right)}^{-1}\right)}^{T}e,d_1^k=\left(\begin{array}{c}{p}^k\\ 0\end{array}\right),d^k=M_2\left(U^k\right)d_0^k+M_3\left(U^k\right)d_1^k$ 获得可行下降方向d^k；

步骤4)由 ${\left(A_k^1\right)}^T\mathrm{\&rho;}=-{\left|\right|d_0^k\left|\right|}^{\mathrm{\&tau;}}e-{\hat{g}}_{L_k}(x^k+d^k)$ 和 ${\left(A_k\right)}^T{\hat{d}}^k={\left(A_k^1\right)}^T{\mathrm{\&rho;}}^k+{\left(A_k^2\right)}^{T}0={\left(A_k^1\right)}^T{\mathrm{\&rho;}}^k$ 获得一个高阶修正方向

步骤5)计算最大值使其满足以下两条件：

$\begin{array}{c}f(U^k+t_k{d}^k+t_k^2{\hat{d}}^k)\&le;f\left(U^k\right)+{\mathrm{\&alpha;t}}_k\&dtri;f\left(U^k\right)d^k,\\ g_j(U^k+t_k{d}^k+t_k^2{\hat{d}}^k)\&le;0,j\&Element;I\end{array};$

步骤6)迭代:

λ^k+1＝λ^k+t_kΔλ^k；

${\mathrm{\&mu;}}^{k+1}\left(i\right)=\min \left\{\max \right\{{\mathrm{\&lambda;}}^{k+1}\left(i\right),{\mathrm{\&lambda;}}^{k0}\left(i\right),\left|\right|d^k\left|\right|\},\stackrel{\&OverBar;}{\mathrm{\&mu;}}\},i\&Element;L_k,$

μ^k+1(i)＝μ^k(i)/2，i∈I/L_k

若积极估计集L_k≠φ，则否则U^k+1＝U^k+t_kd^k；

步骤7)根据方法三更新对称正定矩阵H_k获得H_k+1,k＝k+1,转入步骤2)，直至迭代停止。