CN106803280B - 一种基于变分框架特征感知的细分曲面重建方法 - Google Patents

Abstract

一种基于变分框架的特征感知的细分曲面重建方法，属于数字几何处理和曲面造型领域。首先，根据模型尖锐角、尖锐边和半尖锐边等特征，对控制网格的每一条边设置一个尖锐值，提出了一种特征感知的细分曲面新规则；然后，根据这个新规则，建立基于变分框架和L1范数的细分曲面重建的细分模型；最后，利用增值拉格拉日方法和粒子群优化法迭代求解细分曲面的控制网格的顶点位置和尖锐值，从而得到重建的保特征的细分曲面。该发明提出的基于变分的特征感知的细分曲面重建方法，与现有方法相比，在保持曲面尖锐特征和半尖锐特征等方面具有明显优势，并且能够处理带不同尺度噪音的输入网格的重建问题，在曲面重建上有广泛的应用前景。

Description

一种基于变分框架特征感知的细分曲面重建方法

技术领域

本发明涉及一种基于变分方法的特征感知的细分曲面重建方法，在重建过程中，可以更好地保持曲面尖锐特征、半尖锐特征和角点等几何特征,尤其涉及一种基于变分框架特征感知的细分曲面重建方法，属于数字几何处理和曲面造型领域。

背景技术

细分曲面在计算机图形学中用于从任意网格生成光滑曲面。相对于有理B样条(NURBS)，它具有适用于任意拓扑结构、数值上稳定、实现简易、局部连续性控制等优点。因此，在动画、娱乐、虚拟现实和工业制造等领域有着广泛的应用，并且细分曲面也被集成到MPEG4标准中。

如何保特征的重建细分曲面是细分曲面重建的重要衡量标准，而曲面的几何特征包括尖锐特征(sharp)、半尖锐特征(semi-sharp)和角点(corner)等。细分规则对最终重建的曲面特征有重要的影响，现有的细分规则生成细分曲面主要有两种：一种是生成全局光滑的曲面；另外一种是生成保持尖锐特征的曲面，如Hoppe等人在《Piecewise SmoothSurface Reconstruction》中提出了产生角点、褶皱(crease)和镖点(dart)的细分规则以及重建方法。但是这类方法主要存在两个问题：首先，这类方法需要先检测出具有尖锐特征的边，最终效果非常依赖于初始特征边的检测；其次，此方法控制顶点的边只有两种状态：光滑和尖锐，不能有效地表示半尖锐特征，具有局限性。

DeRose等人在《Subdivision Surface in Character Animation》文中提出现实生活中大部分物体的特征是半尖锐的，因此提出在控制网格每条边上用一个整数控制保特征细分的次数从而控制曲面的弯曲程度，但是这种方法在逆向工程的曲面重建中，优化算法需要求解一个带整数约束的优化方程，非常不方便。

本发明提出了用一个实数来连续控制曲面特征尖锐程度的方法，对控制网格的每条边定义一个尖锐值(sharpness)，通过变分框架，同时寻找最优的控制顶点位置和尖锐值。该方法提出如何用实数控制曲面特征的细分规则，针对输入模型存在不同类型噪音的问题，提出在变分框架中引进L1范数，使得本发明所提方法对不同尺度的噪音不敏感。

发明内容

本发明的目的是克服现有曲面重建技术中存在的不能有效表现半尖锐特征以及优化算法复杂的问题，提出了一种基于变分框架特征感知的细分曲面重建方法。

本发明的核心思想为：用一个实数来连续控制曲面特征尖锐程度，对控制网格的每条边定义一个尖锐值,进一步提出如何用实数控制曲面特征的细分规则；再通过变分框架，同时寻找最优的控制顶点位置和尖锐值；针对输入模型存在不同类型噪音的问题，在变分框架中引进L1范数；首先，根据曲面尖锐角、尖锐边和半尖锐边为主的特征，对控制网格的每一条边设置一个尖锐值，提出了一种特征感知的Loop细分曲面新规则；再根据这个新规则，建立基于变分框架和L1范数的细分曲面重建模型；最后，利用增值拉格拉日方法(Augmented Lagrangian Method)和粒子群优化法(Particle Swarm Optimization)迭代求解细分曲面的控制网格的新的顶点位置和尖锐值，从而得到重建的保特征细分曲面。

一种基于变分框架特征感知的细分曲面重建方法，具体步骤如下：

步骤一：对输入模型进行预处理，得到初始控制网格；

其中，初始控制网格记为

并记输入模型为S⁰，其顶点集合记为

n_v为输入模型顶点个数；

若输入模型不是三角网格模型，则对先对输入模型进行网格重建，将其转化为三角网格模型；

其中，对输入模型的预处理，具体为：

步骤1.1对输入模型采用文献1中3.2节提出的基于ROF模型的去噪方法进行去噪，输出去噪后模型；

文献1：Xiaoqun Wu,Jianmin Zheng,Yiyu Cai,Chi-Wing Fu,Mesh DenoisingUsing Extended ROF Model with L1Fidelity,Computer Graphics Forum,2015,34(7):35-45；

步骤1.2对经过步骤1.1去噪以后输出的去噪后模型，采用文献2中提出的基于二次误差度量的方法进行简化得到初始控制网格

；

文献2：Michal Garland,Paul S.Heckbert,Surface Simplification usingQuadric Error Metrics,Proceeding SIGGRAPH'97Proceedings of the 24th annualconference on Computer Graphics and interactive techniques Pages 209-216.

其中，初始控制网格的顶点集合记为p＝{p₁,p₂,...,p_m}，其中m为控制网格顶点个数，n_e为控制网格边的条数,每一条边都有一个尖锐值，记为e.s∈[0,1]，表明控制多边形的每条边的尖锐及光滑情况，初始控制网格的边可以根据尖锐值进行分类，遵循如下网格边分类原则为：

1.A当e.s＝0时，表明这条边是光滑的；

1.B当e.s＝1时，表明这条边是尖锐的；

1.C当0＜e.s＜1时，表明这条边介于光滑和尖锐之间；

步骤二：将步骤一输出的初始控制进行分类网格细分得到新的边顶点和顶点；

其中，新的边顶点记为e^v，顶点记为v^new；步骤二，具体为：

步骤2.1采用1.A、1.B及1.C描述的网格边分类原则对初始控制网格的每一条边进行分类，得出所有尖锐值组成的控制细分曲面特征；

步骤2.2通过控制细分曲面特征，采用特征感知Loop细分规则生成曲面，步骤为：

步骤2.21，计算新的边顶点e^v；

具体由下面公式(1)得到：

其中，

和

分别是由光滑细分规则和尖锐细分规则产生的边顶点；

新的边顶点由该边的两个端点v₁和v₂，以及与该边相对的两个顶点v₃和v₄这四个顶点加权平均得到，其中，v₁,v₂,v₃,v₄的权值光滑细分规则为

尖锐细分规则为

步骤2.22计算顶点v^new，具体用v.n表示这个顶点的邻域内边尖锐值非零的条数，并根据v.n的大小分如下几种情况：

3A当v.n＜2，这个顶点被认为是光滑的，则该点按照光滑顶点规则计算，即v^new＝v_smooth；光滑规则是由该顶点以及其相邻顶点的加权平均得到，相邻顶点每点的权为

该点的权为(1-a(n))，其中，

n为该顶点相邻的顶点个数；

3B当v.n＝2，这个顶点被认为是褶皱点，假设与其相邻的两条边的尖锐值为零的边分别为e_j和e_k，则该点的尖锐值

则该褶皱点新的位置可由如下公式(2)得到:

其中,v_smooth和v_crease分别为通过光滑规则和尖锐规则得到的顶点位置，褶皱顶点v_crease由该顶点以及其相邻尖锐值非零的边e_j和e_k的另一侧顶点加权平均得到，权分配为

3C当v.n＞2，这个点被认为是角点，假设用A_s表示邻域内的尖锐值非零的边的集合，则这点的尖锐值

则这个新的角点位置可以由如下公式(3)计算得到：

其中,

是由尖锐规则得到的尖锐点的加权平均，v_corner是角点，角点即为该点的位置不变；

步骤三、建立基于变分框架和L1范数的细分曲面重建的能量函数，计算光滑项和误差项，具体为：

步骤3.1，建立基于变分框架和L1范数的细分曲面重建的能量函数；

初始控制网格的顶点集合为p＝{p₁,p₂,...,p_m}，尖锐值集合记为

保持控制顶点的连接关系不变，通过变分框架，可通过优化如下公式(4)所示的细分曲面重建的能量函数得到最优的控制顶点位置和尖锐值；

其中，E_s(p)为光滑项，用来阻止控制网格不出现自交以及不受噪音影响；E_f(p,h)为误差项；

步骤3.2，计算光滑项；光滑项具体通过如下公式(5)计算：

其中，L(p_i)为该点的离散拉普拉斯(Laplacian)，

表示L(p_i)的L₁范数变量i变化范围为1到m，为简化起见，用L表示拉普拉斯矩阵；

步骤3.3，计算误差项；计算输入模型的点到由控制网格生成的极限细分曲面S的距离；

其中，误差项是指用

逼近S，

和S的差；

由

根据步骤二的特征感知的细分规则细分k次得到S^k；按文献3中2.2节提出的规则，对每个顶点得到一个极限点，然后将S^k的所有顶点拉到对应的极限点，

可以表示为初始点p的一个仿射变换,记为M(h)p；为了使得我们的方法不依赖于步骤一预处理中采用的文献1的去噪算法，我们计算输入模型S₀到

的距离；因此误差项计算如下公式(6)：

E_f(p,h)＝||v-M(h)p||＝||v-Mp|| (6)

其中，

为输入模型的顶点集合，M(h)是依赖于h的矩阵，为简化起见，用M表示该矩阵，||||为两个向量差的L₁范数；

综上，细分曲面重建的能量函数(4)可以写成如下公式(7)

文献3：Hugues Hoppe,Tony DeRose,Tom Duchamp,Michael Halstead,HubertJin,John McDonald,Jean Schweitzer,Werner Stuetzle,Piecewise Smooth SurfaceReconstruction,ACM SIGGRAPH 1994Proceedings,295-302

步骤四、应用增值拉格拉日方法求解步骤三的细分曲面重建的能量函数，具体为：

步骤4.1，将(7)转化为如下公式(8)带约束的优化问题：

则根据增值拉格拉日方法可以将上述约束问题转为求如下公式(9)的泛函鞍点问题：

其中，λ_q与λ_z是拉格拉日乘子，<,>是内积，a_p＞0和a_z＞0。则优化问题转化为如下公式(10)鞍点问题：

步骤4.2，求解鞍点问题(10)；将鞍点问题转化为依次求解4个子问题，然后更新拉格拉日乘子的迭代算法，4个子问题分别如下：

步骤4.2A固定h,q,z，求解p子问题，p子问题可以转化如下公式(11)的二次方程形式：

该问题可以转化为线性方程求解；

步骤4.2B固定p,q,z，求解h子问题，h子问题可以转化为如下公式(12)形式：

该问题本方法采用粒子群优化法求解；

步骤4.2C固定p,h,z，求解q子问题，q子问题可以写成如下公式(13)形式：

对公式(13)进行整理，可以转化为如下公式(14)的优化问题：

公式(14)可以分解并且有如下公式(15)的封闭形式解：

其中，

步骤4.2D固定p,h,q，求解z子问题，z子问题有如下公式(16)形式：

则类似(14)，该问题有如下公式(17)的封闭形式解：

其中

步骤4.3，更新拉格拉日乘子，其中第l+1次的迭代为第l次的关系如下：

步骤4.4，迭代求解；

令初值h⁰＝z⁰＝q⁰＝0；

依次迭代求解方程(11)，(12)，(13)，(16)更新拉格拉日乘子(18)，直到满足终止条件；

其中，终止条件为：假设连续两次迭代，如l,l+1次迭代，控制顶点的距离记为

当ε小于给定的阈值ε₀时，迭代停止；

至此，从步骤一到步骤四，完成了一种基于变分框架特征感知的细分曲面重建算法。

有益效果

本发明提出了一种基于变分框架特征感知的细分曲面重建算法，与现有曲面重建方法相比，具有如下有益效果：

1.本发明所提方法在保持曲面尖锐特征和半尖锐特征等方面，我们的算法具有明显优势；

2.本发明所提方法提出了如何用实数控制曲面特征的细分规则，针对输入模型存在不同类型噪音的问题，提出在变分框架中引进L1范数，使得本发明所提方法对不同尺度的噪音不敏感；

3.本发明所提方法在动画产业，虚拟现实和工业制造等领域具有广泛的应用前景。

附图说明

图1是本发明一种基于变分框架特征感知的细分曲面重建算法的框架图；

图2是本发明一种基于变分框架特征感知的细分曲面重建算法及实施例1中Loop细分曲面光滑规则；

图3是本发明一种基于变分框架特征感知的细分曲面重建算法及实施例1中的Loop细分曲面尖锐规则；

图4是本发明一种基于变分框架特征感知的细分曲面重建算法及实施例1中的算法图。

具体实施方式

下面结合附图及实施例对本发明所提方法进行详细说明。

如图1所示，是本发明一种基于变分框架特征感知的细分曲面重建算法及本实施例1的算法流程图；由图1可以看出，本发明的具体实现步骤如下：

步骤A：对输入模型进行预处理，得到初始控制网格

假设输入模型记为S⁰，其顶点集合记为v＝{v₁,v₂,...,v_n}；

若输入模型不是三角网格，则对先对进行网格重建转化为三角网格模型；一般模型都含有噪音，对模型进行去除噪音处理；

对去噪以后的模型，用二次误差度量的方法进行简化得到初始控制网格

步骤B：对输入的初始控制网格和每条边初始的尖锐值e.s∈[0,1]细分得到新的边顶点和顶点；

步骤B.1，计算边顶点e^v。新的边顶点e^v可以由公式(1)计算得到，即：

其中，

和

分别是有光滑规则和尖锐规则产生的边顶点,如图2,图3边蒙板和褶皱边蒙板所示，每个点上的标记的是该点的权重；

步骤B.2，计算顶点v^new。用v.n表示这个顶点的邻域内边尖锐值非零的条数，则根据v.n的大小分如下几种情况：

·当v.n＜2，这个顶点被认为是光滑的，则该点按照光滑顶点规则计算，即v^new＝v_smooth；

·当v.n＝2，这个顶点被认为是褶皱点，假设与其相邻的两条边的尖锐值为零的边分别为e_j和e_k，则该点的的尖锐值

则该褶皱点新的位置可以由公式

得到其中v_smooth和v_crease分别为通过光滑规则和尖锐规则得到的顶点位置，计算如图2和图3中顶点蒙板和褶皱顶点蒙板所示，其中

n是顶点邻域内相邻点的个数；

·当v.n＞2，这个点被认为是角点，假设用A_s表示邻域内的尖锐值非零的边的集合，则这点的尖锐值

则这个新的角点位置可以由

计算得到其中，

是由尖锐规则得到的尖锐点的加权平均，v_corner是角点，计算如图2角点蒙板所示；

步骤C:根据步骤B进行k次特征感知的细分得到细分曲面，建立基于变分框架和L1范数的细分曲面重建的能量函数；

步骤C.1，建立基于变分框架和L1范数的细分曲面重建的能量函数。初始控制网格的顶点集合记为p＝{p₁,p₂,...,p_m}，尖锐值集合记为

其中m为控制网格顶点个数，n_e为控制网格边的条数；保持控制顶点的连接关系不变，通过变分框架，可以通过优化如下能量函数得到最优的控制顶点的位置和尖锐值；

其中E_s(p)为光滑项，E_f(p,h)为误差项；

步骤C.2，计算光滑项；光滑项按

计算；

L(p_i)为该点的离散拉普拉斯(Laplacian)，

为控制顶点离散拉普拉斯的L₁范数和，为简化起见，用L表示拉普拉斯矩阵；

步骤C.3，计算误差项,计算输入模型的点到由控制网格生成的极限细分曲面S的距离；本发明用

逼近S,其中

由

根据步骤B的特征感知的细分规则细分k次得到S^k对每个顶点得到一个极限点，然后将S^k的所有顶点拉到对应的极限点，

可以表示为初始点p的一个仿射变换,记为M(h)p,为了使得我们的方法不依赖于去噪算法，我们计算输入模型到

的距离,因此误差项由

E_f(p,h)＝||v-M(h)p||＝||v-Mp||计算得到。

其中M(h)是依赖于h的矩阵，为简化起见，用M表示该矩阵；

综上，方程

可以写成

步骤D、应用增值拉格拉日方法求解方程

步骤D.1，将方程

转化为如下带约束的优化问题：

则根据增值拉格拉日方法可以将上述约束问题转为求如下泛函的鞍点问题：

其中，λ_q与λ_z是拉格拉日乘子，<,>是内积，a_p＞0和a_z＞0,则优化问题转化为如下鞍点问题：

步骤D.2，求解优化问题

将其转化为依次求解4个子问题，然后更新拉格拉日乘子的迭代算法，4个子问题分别如下：

·固定h,q,z，求解p子问题，p子问题可以转化如下二次方程形式：

该问题可以转化为线性方程求解；

·固定p,q,z，求解h子问题，h子问题可以转化为如下形式：

该问题本方法采用粒子群优化法求解；

·固定p,h,z，求解q子问题，q子问题可以写成如下形式：

对

进行整理，可以转化为如下优化问题：

再将

分解为封闭形式解：

其中，

·固定p,h,q，求解

形式的z子问题，此问题有如下封闭形式解：

其中，

步骤D.3，更新拉格拉日乘子，其中第l+1次的迭代为第l次的关系如下：

步骤D.4，迭代求解；

令初值h⁰＝z⁰＝q⁰＝0；

依次迭代求解方程

和

更新拉格拉日乘子

直到满足终止条件，见图4算法图。

以上所述为本发明的较佳实施例而已，本发明不应该局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。

Claims (5)

1.一种基于变分框架特征感知的细分曲面重建方法，其特征在于：具体步骤如下：

步骤一：对输入模型进行预处理，得到初始控制网格；

其中，初始控制网格记为

并记输入模型为S⁰，其顶点集合记为

n_v为输入模型顶点个数；

若输入模型不是三角网格模型，则对先对输入模型进行网格重建，将其转化为三角网格模型；

其中，对输入模型的预处理，具体为：

步骤1.1对输入模型采用文献1中3.2节提出的基于ROF模型的去噪方法进行去噪，输出去噪后模型；

文献1：Xiaoqun Wu,Jianmin Zheng,Yiyu Cai,Chi-Wing Fu,Mesh Denoising UsingExtended ROF Model with L1 Fidelity,Computer Graphics Forum,2015,34(7):35-45；

步骤1.2对经过步骤1.1去噪以后输出的去噪后模型，采用文献2中提出的基于二次误差度量的方法进行简化得到初始控制网格

文献2：Michal Garland,Paul S.Heckbert,Surface Simplification using QuadricError Metrics,Proceeding SIGGRAPH'97Proceedings of the 24th annual conferenceon Computer Graphics and interactive techniques Pages 209-216.

其中，初始控制网格的顶点集合记为p＝{p₁,p₂,...,p_m}，其中m为控制网格顶点个数，n_e为控制网格边的条数,每一条边都有一个尖锐值，记为e.s∈[0,1]，表明控制多边形的每条边的尖锐及光滑情况，初始控制网格的边根据尖锐值进行分类，遵循如下网格边分类原则为：

1.A当e.s＝0时，表明这条边是光滑的；

1.B当e.s＝1时，表明这条边是尖锐的；

1.C当0＜e.s＜1时，表明这条边介于光滑和尖锐之间；

步骤二：将步骤一输出的初始控制进行分类网格细分得到新的边顶点和顶点；

步骤三、建立基于变分框架和L1范数的细分曲面重建的能量函数，计算光滑项和误差项；

步骤四、应用增值拉格拉日方法求解步骤三的细分曲面重建的能量函数。

2.如权利要求1所述的一种基于变分框架特征感知的细分曲面重建方法，其特征在于：步骤二中，新的边顶点记为e^v，顶点记为v^new；步骤二，具体为：

步骤2.1采用1.A、1.B及1.C描述的网格边分类原则对初始控制网格的每一条边进行分类，得出所有尖锐值组成的控制细分曲面特征；

步骤2.2通过控制细分曲面特征，采用特征感知Loop细分规则生成曲面。

3.如权利要求2所述的一种基于变分框架特征感知的细分曲面重建方法，其特征在于：步骤2.2，具体为：

步骤2.21，计算新的边顶点e^v；

具体由下面公式(1)得到：

其中，

和

分别是由光滑细分规则和尖锐细分规则产生的边顶点；

新的边顶点由该边的两个端点v₁和v₂，以及与该边相对的两个顶点v₃和v₄这四个顶点加权平均得到，其中，v₁,v₂,v₃,v₄的权值光滑细分规则为

尖锐细分规则为

步骤2.22计算顶点v^new，具体用v.n表示这个顶点的邻域内边尖锐值非零的条数，并根据v.n的大小分如下几种情况：

3A当v.n＜2，这个顶点被认为是光滑的，则该点按照光滑顶点规则计算，即v^new＝v_smooth；光滑规则是由该顶点以及其相邻顶点的加权平均得到，相邻顶点每点的权为

该点的权为(1-a(n))，其中，

n为该顶点相邻的顶点个数；

3B当v.n＝2，这个顶点被认为是褶皱点，与其相邻的两条边的尖锐值为零的边分别为e_j和e_k，则该点的尖锐值

则该褶皱点新的位置由如下公式(2)得到:

其中,v_smooth和v_crease分别为通过光滑规则和尖锐规则得到的顶点位置，褶皱顶点v_crease由该顶点以及其相邻尖锐值非零的边e_j和e_k的另一侧顶点加权平均得到，权分配为

3C当v.n＞2，这个点被认为是角点，用A_s表示邻域内的尖锐值非零的边的集合，则这点的尖锐值

则这个新的角点位置由如下公式(3)计算得到：

其中,

是由尖锐规则得到的尖锐点的加权平均，v_corner是角点，角点即为该点的位置不变。

4.如权利要求1所述的一种基于变分框架特征感知的细分曲面重建方法，其特征在于：步骤三，具体为：

步骤3.1，建立基于变分框架和L1范数的细分曲面重建的能量函数；

初始控制网格的顶点集合为p＝{p₁,p₂,...,p_m}，尖锐值集合记为

保持控制顶点的连接关系不变，通过变分框架，通过优化如下公式(4)所示的细分曲面重建的能量函数得到最优的控制顶点位置和尖锐值；

其中，E_s(p)为光滑项，用来阻止控制网格不出现自交以及不受噪音影响；E_f(p,h)为误差项；

步骤3.2，计算光滑项；光滑项具体通过如下公式(5)计算：

其中，L(p_i)为该点的离散拉普拉斯Laplacian，

表示L(p_i)的L₁范数变量i变化范围为1到m，为简化起见，用L表示拉普拉斯矩阵；

步骤3.3，计算误差项；计算输入模型的点到由控制网格生成的极限细分曲面S的距离；

其中，误差项是指用

逼近S，

和S的差；

由

根据步骤二的初始控制进行分类网格细分，细分k次得到S^k；按文献3中2.2节提出的规则，对每个顶点得到一个极限点，然后将S^k的所有顶点拉到对应的极限点，

表示为初始点p的一个仿射变换,记为M(h)p；计算输入模型S⁰到

的距离；因此误差项计算如下公式(6)：

E_f(p,h)＝||v-M(h)p||＝||v-Mp|| (6)

其中，

为输入模型的顶点集合，M(h)是依赖于h的矩阵，为简化起见，用M表示该矩阵，||||为两个向量差的L₁范数；

综上，将细分曲面重建的能量函数(4)写成如下公式(7):

文献3：Hugues Hoppe,Tony DeRose,Tom Duchamp,Michael Halstead,Hubert Jin,John McDonald,Jean Schweitzer,Werner Stuetzle,Piecewise Smooth SurfaceReconstruction,ACM SIGGRAPH 1994 Proceedings,295-302。

5.如权利要求4所述的一种基于变分框架特征感知的细分曲面重建方法，其特征在于：步骤四，具体为：

步骤4.1，将(7)转化为如下公式(8)带约束的优化问题：

则根据增值拉格拉日方法将上述约束问题转为求公式(9)的泛函鞍点问题：

其中，λ_q与λ_z是拉格拉日乘子，<,>是内积，a_p＞0和a_z＞0，则优化问题转化为公式(10)鞍点问题：

步骤4.2，求解鞍点问题(10)；将鞍点问题转化为依次求解4个子问题，然后更新拉格拉日乘子的迭代算法，4个子问题分别如下：

步骤4.2A固定h,q,z，求解p子问题，p子问题转化如下公式(11)的二次方程形式：

该问题转化为线性方程求解；

步骤4.2B固定p,q,z，求解h子问题，h子问题转化为如下公式(12)形式：

该问题本方法采用粒子群优化法求解；

步骤4.2C固定p,h,z，求解q子问题，q子问题写成如下公式(13)形式：

对公式(13)进行整理，转化为如下公式(14)的优化问题：

公式(14)分解并且有如下公式(15)的封闭形式解：

其中，

步骤4.2D固定p,h,q，求解z子问题，z子问题有如下公式(16)形式：

则类似(14)，该问题有如下公式(17)的封闭形式解：

其中

步骤4.3，更新拉格拉日乘子，其中第l+1次的迭代为第l次的关系如下公式(18)：

步骤4.4，迭代求解；

令初值h⁰＝z⁰＝q⁰＝0；

依次迭代求解方程(11)，(12)，(13)，(16)更新拉格拉日乘子(18)，直到满足终止条件；

其中，终止条件为：假设连续两次l,l+1次迭代，控制顶点的距离记为

当ε小于给定的阈值ε₀时，迭代停止。

Images