CN107229595B - 一种soa架构的电网实时抗差状态估计方法 - Google Patents

Abstract

一种SOA架构的电网实时抗差状态估计系统及方法。解决电力导读系统中数据传输会渗入噪音，影响系统正确计算的问题。系统，包括相连的用户端、服务器端，用户端与电网系统相连接。服务端包括处理单元、服务接口单元，用户端与服务接口单元通过网络相连接。方法对采集的量测和拓扑信息采用最小绝对值法进行状态估计计算，有效降低了电力调度自动化系统中原始数据所含的噪声及剔除其中的不良数据，保证了电力调度系统软件计算的准确性。

Description

一种SOA架构的电网实时抗差状态估计方法

技术领域

本发明涉及一种电力调度自动化技术领域，尤其是涉及一种SOA 架构的电网实时抗差状态估计方法。

背景技术

电力系统中的运行数据从变电站RTU采集到经过网络传输到电力调度自动化系统中时，往往会受到随机干扰而渗入误差，甚至会产生严重不良数据，这些数据也称为生数据。生数据往往会影响调度系统高级应用软件的准确计算。因此，降低电力调度自动化系统中生数据所含的噪声及剔除其中的不良数据对于调度系统的正常运行至关重要。

此外，当电网中出现因不同电源供电而造成电网分为多区域独立运行的情况时，在调度自动化系统中，通常状态估计方法需要对不同的区域分别进行独立计算，将就会增加计算流程和时间。因此，发明了一种将所有独立区域的状态估计计算过程合并为一个计算过程的计算方法，以达到简化计算流程、减少计算时间的效果。

发明内容

本发明主要是解决电力调度自动化系统中的生数据传输含有噪声或不良数据，影响系统高级应用软件准确计算的问题，提供了一种 SOA架构的电网实时抗差状态估计方法。

本发明的上述技术问题主要是通过下述技术方案得以解决的：一种SOA架构的电网实时抗差状态估计系统，包括

部署在云计算平台上对状态估计系统进行设置和结果显示的用户端，用户端与电网系统相连接，接收来自电力系统的量测值和网架拓扑信息，将这两类数据编辑为标准格式的数据文件，并将数据文件传输给服务端进行状态估计计算；其中量测值包括节点注入功率量测值、支路功率量测和电压幅值量测；网架拓扑信息包括系统节点数、平衡节点、线路支路参数和变压器支路参数。

部署在云计算平台上用于实现状态估计逻辑计算的服务端，服务端包括根据电网量测值和网架拓扑信息进行状态估计计算的处理单元，以及供用户端接收数据的服务接口单元，用户端与服务接口单元通过网络相连接。本发明中用户端和服务端均部署在云计算平台，与传统的双机热备用模式相比，不仅节约了成本，更重要的是提高了应用本身的运行可靠性和应对各种灾害的能力，具有自动灾备功能。状态估计系统采用SOA构架，以面向服务的形式搭建，使得应用具有松散耦合、跨平台和可复用等优点，与电网调度自动化系统的结合具有高度灵活性。并采用B/S模式，将计算的结果以网页形式显示，用户端只需浏览器便可进行访问。状态估计应用的界面中的主网图采用 SVG可伸缩矢量图形技术，主网图显示了实时遥测数据和实时状态估计结果，以及显示了状态估计计算过程中的收敛曲线及各量测值对应的计算值。服务端采用IIS托管。本发明服务端通过状态估计计算，有效降低了电力调度自动化系统中生数据所含的噪声及剔除其中的不良数据。

一种SOA架构的电网实时抗差状态估计方法，包括以下步骤：

S1.用户端获取电网量测值和网架拓扑信息；

S2.根据用户对量测值的选择调整作为状态估计计算的原始数据；

S3.根据原始数据分析电网在电气上是否为单一区域运行情况，若是，则直接将原始数据处理生成一个区域的标准格式数据文件，若不是，则按照多区域处理办法生成合并计算的标准格式数据文件；

S4.将标准格式的数据文件发送给服务端；

S5.服务端读取数据文件中的数据进行状态估计计算，计算结果供用户端进行调用。状态估计应用服务端内嵌了五种特点各异的状态估计方法，用户可根据实际电网情况选择一种运行，本发明优先采用最小绝对值法。

作为一种优选方案，步骤S1中获取的电网量测值包括节点注入功率量测值、节点电压量测值、支路功率量测值、零注入功率量测值；网架拓扑信息包括电网节点数、平衡节点号、节点连接关系、线路和变压器支路阻抗参数、接地支路参数。

作为一种优选方案，步骤S2的具体过程包括：

S21.将电网中所有可以选择的量测数据点在用户端上显示并供用户选择；

S22.根据用户选择的量测数据点并结合默认选择的量测数据点，与上一次进行状态估计计算的量测数据点进行比较；

S23.将上一次进行状态估计计算的量测数据点作为原始数据，若比较为增加了量测数据点，将增加的数据点添加到原始数据中，若比较为减少了量测数据点，将减少的数据点从原始数据中删除，若比较为没有改变，则保持原始数据量测点不变。

作为一种优选方案，步骤3中计算标准格式数据文件的具体过程包括：

S31.根据网架拓扑连接关系识别电网在电气关系上是一个电气区域，还是相互独立的多个电气区域，即是否为多区域独立运行情况；

S32.若是一个区域，则按常规方法将原始数据处理为标准格式数据文件，进入步骤S4；若为多个区域，则将每个区域的原始数据分别处理为标准格式数据文件；

S33.以节点最多的电气区域为基础，保持该区域所有节点的编号不变，其他区域的节点号依次按顺序递增，并重新按照标准格式数据文件进行处理，直到将所有区域的标准格式数据合并为一个多平衡节点的标准格式数据文件；假如整个电网有100个节点，分成3个独立的区域运行，区域1、2和3的节点数分别为50、30和20。则将节点数为50的区域1设为基础区域，将其所有节点分别编号为1～50，将区域2的节点编号为51～80，将区域3的节点编号为81～100。若电网分为更多的区域，以此类推进行处理。

S34.合并后，按常规方法将原始数据处理为标准格式数据文件。本方案采用对电网中所有独立运行的区域合并计算的方法，计算结果与独立计算结果一致，避免了计算重复性，有效简化了计算流程，提高了计算效率。

该标准格式数据文件的具体数据结构为：

1)节点数，支路数，基准功率，迭代次数，中心参数；

2)收敛精度；

3)序号，平衡节点1，平衡节点2，…；

4)线路参数(序号，节点i，节点j，r，x，b/2)；

5)变压器参数(序号，节点i，节点j，r，x，k0)；

6)接地支路参数(节点i，gi，bi)；

7)零注入量测(序号，节点i)；

8)节点零注入量测(序号，节点i，Pgi，Qgi，Pdi，Qdi)；

9)支路量测(序号，节点i，节点j，Pij，Qij)；

10)电压量测(序号，节点i，Vi)；

11)电压运行限制(序号，节点i，下限，上限)；

12)有功出力限制(序号，节点i，下限，上限)；

13)无功出力限制(序号，节点i，下限，上限)；

如果为电网多区域运行。则3)处有多个平衡节点，若只有一个网架，则只有平衡节点1；r表示电阻，x表示电抗，b表示电纳，k0 表示变压器当前变比。

作为一种优选方案，步骤5中状态估计计算的具体过程包括：

S51.读取数据文件，根据状态估计的非线性量测方程

z＝g(x)+ε

采用最小绝对值法进行计算，最小绝对值模型为

其中z为m维量测向量，即三类量测值中用于计算的m各量测值组成的向量，z＝[P_i,Q_i；V_i；P_ij,Q_ij]，P_i、Q_i分别为节点注入有功率和无功功率量测值，V_i为节点电压幅值量测值，P_ij、Q_ij为支路功率量测值；x为n维状态向量，即x＝[V,δ]，V、δ分别为电网节点电压幅值和相角，共n个；g(x)为m维量测函数向量，为量测值z中m 个量测值对应的m个函数表达式向量，即为估计值，

为节点注入功率估计值，

为节点电压估计值，

为支路功率估计值；ε为m维量测残差，h(x)为节点零节点注入功率量测值，w为m维权重向量；T为矩阵运算符号，表示对向量的转置；

S52.初始化参数，设定迭代次数k＝0，k_max＝20中心参数σ＝0.1， l＝u＝λ＝[1,…,1]^T∈R^m，y＝0，α＝β＝w，收敛精度e＝10^-5，状态变量x采用平启动，即V＝1，δ＝0，其中l,u为松弛变量， (l,u)∈R^m≥0，y、λ、α、β为拉格朗日乘子向量， y∈R^c，λ、α、β∈R^m，R^c、R^m分别为c维和m维实数向量；

S53.判断是否k<k_max，若是，进入下一步骤，若否，最小绝对值模型不收敛，计算结束；

S54.计算互补间隙Gap＝α^Tl+β^Tu，若Gap<e，则最小绝对值模型收敛，输出结果x，若否继续下一步骤；

S55.计算扰动因子

S56.根据最小绝对值模型获得修正方程，求解修正方程得到修正量Δx，Δλ，Δy，Δα，Δβ，Δl，Δu；

S57.求出原始对偶步长

i＝1,…,m；

S58.修正原始、对偶变量

将k＝k+1，返回步骤S53 。

作为一种优选方案，步骤S56 修正量获取的过程包括：

S561.等价转化最小绝对值模型，获得最小值等价模型

引入松弛变量(l,u)∈R^m≥0将不等式转化为等式，并进行等价转换，得到

S562.构建拉格朗日函数：

L＝w^T(l+u)-y^Th(x)-λ^T(z-g(x)+l-u)-α^Tl-β^Tu (3)

求取拉格朗日函数极限，即为所求状态估计结果，导出下一阶最优性条件：

L_y≡h(x)＝0 (5)

L_λ≡z-g(x)+l-u＝0 (6)

L_l≡w-λ-α＝0 (7)

L_u≡w+λ-β＝0 (8)

L_α≡diag(α)diag(l)e＝0 (9)

L_β≡diag(β)diag(u)e＝0 (10)

其中，e＝[1,…,1]^T∈R^m，引入扰动参数

对互补项(9)、 (10)实施扰动松弛得到：

S563.对以上非线性方程采用牛顿法求解，线性化得到：

-Δλ-Δα＝-L_l (16)

Δλ-Δβ＝-L_u (17)

S564.变量α，β，λ选取合适初始值，使得L_l＝L_u＝0，由公式 (16)、(17)、(18)、(19)求得以下关系式：

求得Δα，Δβ，Δl，Δu；

S565.将Δl、Δu代入公式(15)，得到：

公式中S＝diag(α)^-1diag(l)+diag(β)^-1 diag(u)；

由公式(13)、(14)、(21)得到修正方程：

公式中

求解修正方程得到Δx，Δλ，Δy。

因此，本发明的优点是：

1.通过状态估计计算，有效降低了电力调度自动化系统中生数据所含的噪声及剔除其中的不良数据，保证了电力调度系统软件计算的准确性。

2.采用对电网中所有独立运行的区域合并计算的方法，计算结果与独立计算结果一致，避免了计算重复性，有效简化了计算流程，提高了计算效率。

3.采用SOA架构搭建状态估计系统，有效解决了传统软件架构紧耦合的问题，使得状态估计系统具备松散耦合特性，灵活性高，能跨平台且功能复用，提高了系统的可拓展性。

4.采用B/S模式，克服了C/S模式对单机内存的高要求和固定地点办公的缺点，使得用户使用轻型设备通过互联网便可访问状态估计系统，且具备移动办公功能。

5.SOA架构的状态估计系统部署于云计算平台，克服了传统双机热备用抗灾能力不足的缺点，具有异地灾备能力，能抵抗例如台风、地震等自然灾害；还充分利用了云计算的弹性计算技术，能快速地扩充计算能力。

附图说明

附图1是本发明中状态估计计算的流程示意图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

本实施例一种SOA架构的电网实时抗差状态估计系统，包括

部署在云计算平台上对状态估计系统进行设置和结果显示的用户端，用户端与电网系统相连接，接收来自电力系统的量测值和网架拓扑信息，

部署在云计算平台上用于实现状态估计逻辑计算的服务端，服务端包括根据电网量测值和网架拓扑信息进行状态估计计算的处理单元，以及供用户端接收计算结果的服务接口单元，用户端与服务接口单元通过网络相连接。

一种SOA架构的电网实时抗差状态估计方法，包括以下步骤：

S1.用户端获取电网量测值和网架拓扑信息；

电网量测值包括节点注入功率量测值、节点电压量测值、支路功率量测值、零注入功率量测值；网架拓扑信息包括电网节点数、平衡节点号、节点连接关系、线路和变压器支路阻抗参数、接地支路参数。

S2.根据用户对量测值的选择调整作为状态估计计算的原始数据；具体过程为：

S21.将电网中所有可以选择的量测数据点在用户端上显示并供用户选择；

S22.根据用户选择的量测数据点并结合默认选择的节点量测数据点，与上一次进行状态估计计算的量测数据点进行比较；

S23.将上一次进行状态估计计算的量测数据点作为原始数据，若比较为增加了量测数据点，将增加的数据点添加到原始数据中，若比较为减少了量测数据点，将减少的数据点从原始数据中删除，若比较为没有改变，则保持原始数据量测点不变。

S3.根据原始数据分析电网在电气上是否为单一区域运行情况，若是，则直接将原始数据处理生成一个区域的标准格式数据文件，若不是，则按照多区域处理办法生成合并计算的标准格式数据文件；具体过程为：

S31.根据网架拓扑连接关系识别电网在电气关系上是一个电气区域，还是相互独立的多个电气区域，即是否为多区域独立运行情况；

S32.若是一个区域，则按常规方法将原始数据处理为标准格式数据文件，进入步骤S4；若为多个区域，则将每个区域的原始数据分别处理为标准格式数据文件；

S33.以节点最多的电气区域为基础，保持该区域所有节点的编号不变，其他区域的节点号依次按顺序递增，并重新按照标准格式数据文件进行处理，直到将所有区域的标准格式数据合并为一个多平衡节点的标准格式数据文件；

S34.合并后，按常规方法将原始数据处理为标准格式数据文件。

本实施例以IEEE4和IEEE14节点系统为例。IEEE4和IEEE14 节点系统为独立区域，仿真多区域合并计算方法，计算如下：

设IEEE14节点系统为区域1，标准格式数据文件为：

设IEEE4节点系统为区域2，标准格式数据文件为：

区域1和区域2合并后，标准格式数据文件为：

，则IEEE14和IEEE4节点独立计算和合并计算迭代次数及分别进行 10次的平均耗时如下表所示：

。从表中可知区域1和区域2合并计算后，比独立计算具有更高的效率。合并计算使得多个区域的网架能一次性完成计算，减少了基础计算量。对于规模较小的系统，合并计算能有效减少计算实际，具有一定的工程应用价值。

S4.将标准格式的数据文件发送给服务端；

S5.服务端读取数据文件中的数据进行状态估计计算，计算结果供用户端进行调用。

如图1所示，状态估计计算过程包括：

S51.读取数据文件，根据状态估计的非线性量测方程

z＝g(x)+ε

采用最小绝对值法进行计算，最小绝对值模型为

其中z为m维量测向量，即三类量测值中用于计算的m各量测值组成的向量，z＝[P_i,Q_i；V_i；P_ij,Q_ij]，P_i、Q_i分别为节点注入有功率和无功功率量测值，V_i为节点电压幅值量测值，P_ij、Q_ij为支路功率量测值；x为n维状态向量，即x＝[V,δ]，V、δ分别为电网节点电压幅值和相角，共n个；g(x)为m维量测函数向量，为量测值z中m 个量测值对应的m个函数表达式向量，即为估计值，

为节点注入功率估计值，

为节点电压估计值，

为支路功率估计值；ε为m维量测残差，h(x)为节点零节点注入功率量测值，w为m维权重向量；T为矩阵运算符号，表示对向量的转置；

S52.初始化参数，设定迭代次数k＝0，k_max＝20中心参数σ＝0.1， l＝u＝λ＝[1,…,1]^T∈R^m，y＝0，α＝β＝w，收敛精度e＝10^-5，状态变量x采用平启动，即V＝1，δ＝0，其中l,u为松弛变量， (l,u)∈R^m≥0，y、λ、α、β为拉格朗日乘子向量， y∈R^c，λ、α、β∈R^m，R^c、R^m分别为c维和m维实数向量；

S53.判断是否k<k_max，若是，进入下一步骤，若否，最小绝对值模型不收敛，计算结束；

S54.计算互补间隙Gap＝α^Tl+β^Tu，若Gap<e，则最小绝对值模型收敛，输出结果x，若否继续下一步骤；

S55.计算扰动因子

S56.根据最小绝对值模型获得修正方程，求解修正方程得到修正量Δx，Δλ，Δy，Δα，Δβ，Δl，Δu；其过程包括：

S561.等价转化最小绝对值模型，获得最小值等价模型

引入松弛变量(l,u)∈R^m≥0将不等式转化为等式，并进行等价转换，得到

S562.构建拉格朗日函数：

L＝w^T(l+u)-y^Th(x)-λ^T(z-g(x)+l-u)-α^Tl-β^Tu (3)

求取拉格朗日函数极限，导出下一阶最优性条件：

L_y≡h(x)＝0 (5)

L_λ≡z-g(x)+l-u＝0 (6)

L_l≡w-λ-α＝0 (7)

L_u≡w+λ-β＝0 (8)

L_α≡diag(α)diag(l)e＝0 (9)

L_β≡diag(β)diag(u)e＝0 (10)

其中，e＝[1,…,1]^T∈R^m，引入扰动参数

对互补项(9)、 (10)实施扰动松弛得到：

S563.对以上非线性方程采用牛顿法求解，线性化得到：

-Δλ-Δα＝-L_l (16)

Δλ-Δβ＝-L_u (17)

S564.变量α，β，λ选取合适初始值，使得L_l＝L_u＝0，由公式(16)、 (17)、(18)、(19)求得以下关系式：

求得Δα，Δβ，Δl，Δu；

S565.将Δl、Δu代入公式(15)，得到：

公式中S＝diag(α)^-1diag(l)+diag(β)^-1 diag(u)；

由公式(13)、(14)、(21)得到修正方程：

公式中

求解修正方程得到Δx，Δλ，Δy。

S57.求出原始对偶步长

i＝1,…,m；

S58.修正原始、对偶变量

将k＝k+1，返回步骤S53 。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

尽管本文较多地使用了用户端，服务端，处理单元、网页单元等术语，但并不排除使用其它术语的可能性。使用这些术语仅仅是为了更方便地描述和解释本发明的本质；把它们解释成任何一种附加的限制都是与本发明精神相违背的。

Claims (5)

1.一种SOA架构的电网实时抗差状态估计方法，采用一种SOA架构的电网实时抗差状态估计系统，系统包括部署在云计算平台上对状态估计系统进行设置和结果显示的用户端，用户端与电网系统相连接，接收来自电力系统的量测值和网架拓扑信息，部署在云计算平台上用于实现状态估计逻辑计算的服务端，服务端包括根据电网量测值和网架拓扑信息进行状态估计计算的处理单元，以及供用户端接收数据的服务接口单元，用户端与服务接口通过网络相连接，其特征是包括以下步骤：

S1.用户端获取电网量测值和网架拓扑信息；

S2.根据用户对量测值的选择调整作为状态估计计算的原始数据；

S3.根据原始数据分析电网在电气上是否为单一区域运行情况，若是，则直接将原始数据处理生成一个区域的标准格式数据文件，若不是，则按照多区域处理办法生成合并计算的标准格式数据文件；计算标准格式数据文件的具体过程包括：

S31.根据网架拓扑连接关系识别电网在电气关系上是一个电气区域，还是相互独立的多个电气区域，即是否为多区域独立运行情况；

S32.若是一个区域，则按常规方法将原始数据处理为标准格式数据文件，进入步骤S4；若为多个区域，则将每个区域的原始数据分别处理为标准格式数据文件；

S33.以节点最多的电气区域为基础，保持该区域所有节点的编号不变，其他区域的节点号依次按顺序递增，并重新按照标准格式数据文件进行处理，直到将所有区域的标准格式数据合并为一个多平衡节点的标准格式数据文件；

S34.合并后，按常规方法将原始数据处理为标准格式数据文件；

S4.将标准格式的数据文件发送给服务端；

S5.服务端读取数据文件中的数据进行状态估计计算，计算结果供用户端进行调用。

2.根据权利要求1所述的一种SOA架构的电网实时抗差状态估计方法，其特征是步骤S1中获取的电网量测值包括节点注入功率量测值、节点电压量测值、支路功率量测值、零注入功率量测值；网架拓扑信息包括电网节点数、平衡节点号、节点连接关系、线路和变压器支路阻抗参数、接地支路参数。

3.根据权利要求1所述的一种SOA架构的电网实时抗差状态估计方法，其特征是步骤S2的具体过程包括：

S21.将电网中所有可以选择的量测数据点在用户端上显示并供用户选择；

S22.根据用户选择的量测数据点并结合默认选择的节点量测数据点，与上一次进行状态估计计算的量测数据点进行比较；

S23.将上一次进行状态估计计算的量测数据点作为原始数据，若比较为增加了量测数据点，将增加的数据点添加到原始数据中，若比较为减少了量测数据点，将减少的数据点从原始数据中删除，若比较为没有改变，则保持原始数据量测点不变。

4.根据权利要求1或2所述的一种SOA架构的电网实时抗差状态估计方法，其特征是步骤5中状态估计计算的具体过程包括：

S51.读取数据文件，根据状态估计的非线性量测方程

z＝g(x)+ε

采用最小绝对值法进行计算，最小绝对值模型为

其中z为m维量测向量，即三类量测值中用于计算的m各量测值组成的向量，z＝[P_i,Q_i；V_i；P_ij,Q_ij]，P_i、Q_i分别为节点注入有功率和无功功率量测值，V_i为节点电压幅值量测值，P_ij、Q_ij为支路功率量测值；x为n维状态向量，即x＝[V,δ]，V、δ分别为电网节点电压幅值和相角，共n个；g(x)为m维量测函数向量，为量测值z中m个量测值对应的m个函数表达式向量，即为估计值，

为节点注入功率估计值，

为节点电压估计值，

为支路功率估计值；ε为m维量测残差，h(x)为节点零节点注入功率量测值，w为m维权重向量；T为矩阵运算符号，表示对向量的转置；

S52.初始化参数，设定迭代次数k＝0，k_max＝20中心参数σ＝0.1，l＝u＝λ＝[1,…,1]^T∈R^m，y＝0，α＝β＝w，收敛精度e＝10^-5，状态变量x采用平启动，即V＝1，δ＝0，其中l,u为松弛变量，(l,u)∈R^m≥0，y、λ、α、β为拉格朗日乘子向量，y∈R^c，λ、α、β∈R^m，R^c、R^m分别为c维和m维实数向量；

S53.判断是否k<k_max，若是，进入下一步骤，若否，最小绝对值模型不收敛，计算结束；

S54.计算互补间隙Gap＝α^Tl+β^Tu，若Gap<e，则最小绝对值模型收敛，输出结果x，若否继续下一步骤；

S55.计算扰动因子

S56.根据最小绝对值模型获得修正方程，求解修正方程得到修正量Δx，Δλ，Δy，Δα，Δβ，Δl，Δu；

S57.求出原始对偶步长

i＝1,…,m；

S58.修正原始、对偶变量

将k＝k+1，返回步骤S53 。

5.根据权利要求1所述的一种SOA架构的电网实时抗差状态估计方法，其特征是步骤S56 修正量获取的过程包括：

S561.等价转化最小绝对值模型，获得最小值等价模型

引入松弛变量(l,u)∈R^m≥0将不等式转化为等式，并进行等价转换，得到

S562.构建拉格朗日函数：

L＝w^T(l+u)-y^Th(x)-λ^T(z-g(x)+l-u)-α^Tl-β^Tu (3)

求取拉格朗日函数极限，导出下一阶最优性条件：

L_y≡h(x)＝0 (5)

L_λ≡z-g(x)+l-u＝0 (6)

L_l≡w-λ-α＝0 (7)

L_u≡w+λ-β＝0 (8)

L_α≡diag(α)diag(l)e＝0 (9)

L_β≡diag(β)diag(u)e＝0 (10)

其中，e＝[1,…,1]^T∈R^m，引入扰动参数μ＞0，对互补项(9)、(10)实施扰动松弛得到：

S563.对以上非线性方程采用牛顿法求解，线性化得到：

-Δλ-Δα＝-L_l (16)

Δλ-Δβ＝-L_u (17)

S564.变量α，β，λ选取合适初始值，使得L_l＝L_u＝0，由公式(16)、(17)、(18)、(19)求得以下关系式：

求得Δα，Δβ，Δl，Δu；

S565.将Δl、Δu代入公式(15)，得到：

公式中S＝diag(α)^-1diag(l)+diag(β)^-1diag(u)；

由公式(13)、(14)、(21)得到修正方程：

公式中

求解修正方程得到Δx，Δλ，Δy。

Images