张瑞友, 王超慧, 陈勇强
东北大学 信息科学与工程学院, 辽宁 沈阳 110819
收稿日期：2020-02-18
基金项目：国家自然科学基金资助项目(71971050，71831006)。
作者简介：张瑞友(1979-)，男，辽宁沈阳人，东北大学教授，博士生导师。

摘要：为了高效求解非线性规划问题, 对一种基于控制思想的新颖方法——李雅普诺夫方法——进行了研究.该方法将约束非线性规划问题转化为一个动态系统, 基于系统的动态特性给出原优化问题的最优解.分别针对单目标和多目标的非线性规划问题, 对算法的收敛性进行了分析, 给出了算法在应用时松弛变量、增益因子等关键参数的取值建议.大量数值算例验证了上述收敛性及参数取值建议的正确性, 表明了该方法在求解非线性规划问题时的巨大潜力和新颖性.
关键词：约束非线性规划多目标优化李雅普诺夫方法动态系统最优化算法
Lyapunov Method for Solving Nonlinear Programming Problems Based on Control Ideas
ZHANG Rui-you, WANG Chao-hui, CHEN Yong-qiang
School of Information Science & Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: ZHANG Rui-you, E-mail: zhangruiyou@ise.neu.edu.cn.

Abstract: In order to solve nonlinear programming problems efficiently, a novel optimization method named Lyapunov theory-based method (for short, Lyapunov method) based on control ideas is studied. This method transforms a constrained nonlinear programming problem into a dynamic system and presents optimal solution of the original optimization problem according to the dynamic characteristics of the system. Regarding to the single-objective and multi-objective nonlinear programming problems, the convergence of the algorithm is analyzed, and potential values of the key parameters such as the slack variables and gain factors in applications of the algorithm are suggested. A large number of numerical instances verify the aforementioned convergence and the correctness of the proposed parameter values, indicating the great potential and novelty of the method in solving nonlinear programming problems.
Key words: constrained nonlinear programmingmulti-objective optimizationLyapunov methoddynamic systemoptimization algorithm
基于李雅普诺夫理论的方法(简称李雅普诺夫方法)是解决约束非线性规划问题的一种新颖方法.约束非线性规划问题广泛存在于智能制造等领域, 多年来一直是学术界的研究热点与难点之一, 解决这类问题的方法主要有牛顿迭代法等启发式方法^[1]及各种智能优化方法(或称为元启发式算法)^[2].对于单目标非线性规划问题的求解, Hu等^[3]对正交约束问题设计了一种结构拟牛顿算法进行求解.Liu等^[4]对等式约束的最优化问题提出结合线性搜索的序列二次规划方法, 建立全局和局部收敛性理论.Hong等^[5]及Li等^[6]提出用交替方向乘子法求解一类非凸的问题.这些智能优化方法通常计算时间较长且难以保证收敛.
对于多目标非线性规划问题的求解, 进化算法作为一类启发式搜索算法, 已被成功应用于多目标优化领域.首先提出的算法有MOGA和NSGA^[7], 将进化算法与多目标优化问题有机结合.为了提高算法的效率, 又出现了NSGA-II^[8], SPEA2和PESA-II^[9]等算法, 但对于较复杂的多目标优化问题, 这些进化算法仍存在一些缺点, 如难以收敛到理想Pareto前沿、难以获得分布均匀的Pareto最优解、算法停止准则不明确等.
李雅普诺夫方法具较强的处理约束能力和收敛性.李雅普诺夫方法将约束优化问题描述为一个动态系统, 进而基于动态系统的稳定控制状态实现对优化问题的求解.Torelli等^[10]最早于2013年提出李雅普诺夫方法, 将非线性规划问题描述为一个偏微分方程组, 针对无约束的情形分析了该方法的收敛性, 并针对电力系统中的最优潮流问题进行初步的测试.Torelli等^[11]借鉴李雅普诺夫方法将约束转化为微分方程的思想, 提出一种求解微分代数方程的新方法.Torelli等^[12]指出: 牛顿法等传统迭代算法存在因一阶条件的高度非线性而无法收敛等局限, 而李雅普诺夫方法能够克服这些局限性.Torelli等^[13]把二阶导数信息引入李雅普诺夫方法, 进一步提高该方法的鲁棒性.Liang等^[14]证明李雅普诺夫方法中动力系统平衡点的全局收敛性, 并将该方法与牛顿-拉普拉斯算法结合进而求解非线性问题.
综上所述, 李雅普诺夫方法是近年来提出的一种新颖的方法, 在求解约束非线性规划问题方面具有巨大的发展潜力.然而, 针对该算法的研究刚刚开始, 对算法的收敛性和参数设置等研究仍远不完善.本文的主要贡献如下: 1)对李雅普诺夫方法求解单目标约束非线性规划问题和多目标约束非线性规划问题的收敛性进行了理论分析, 得出了若干结论; 2)提出了松弛变量、决策变量增益因子、松弛变量增益因子等关键参数的取值建议, 并给出了数值算例.
1 李雅普诺夫方法的基本思想李雅普诺夫方法求解非线性规划问题的基本过程如下所示^[12].
任意给出下面形式的单目标非线性优化问题(P1):

(1)

(2)

(3)

其中: x为决策变量；f为目标函数; g(x) 为n维等式约束；h(x) 为m维不等式约束, h_min和h_max分别为不等式约束的下界和上界.
引入非负松弛变量q(≥0), 将目标函数(1)的最小化转化为如下偏差函数，

(4)

的零点问题.类似地, 引入非负松弛变量s(≥0)和t(≥0), 问题P1的约束条件(2)和(3)可以转化为

(5)

(6)

(7)

的零点问题.
方程组

(8)

(9)

的解为问题P1的可行解.将方程组(8), (9)进一步转换为标量半正定函数最小化问题:

(10)

其中, E^T (z)表示对E(z)求转置.当E(z^*)=0时, 式(10)有最小值, 即W(z^*)=0, z^*为式(8)的解.
求解式(10)的最小化问题, 也就是求解如下方程：

(11)

由于问题的非线性特性，式(11)难以直接求解，设计1个动态系统, 其中的决策变量和松弛变量与时间相关.在此假设下, 可将式(10)视为李雅普诺夫函数.根据李雅普诺夫定理^[13]可知: 若能使式(10)的时间导数为负定或半负定, 则该动态系统可以稳定收敛于平衡点.
由文献[2]可知, W(z) 对时间的求导结果为

(12)

其中表示对W(z) 求时间导数.为了保证为负定或半负定, 令z的时间导数为

(13)

其中k=(k_x, k_s, k_t, k_q), k_x, k_s, k_t, k_q均为正增益因子, 此时可得

(14)

此时, W (z) 为半正定, 为半负定, 符合李雅普诺夫定理, 系统渐进稳定, 存在平衡点.动态系统平衡时对应的解也就是问题P1的可行解.解的质量与动态系统的平衡过程、平衡位置直接相关.
2 单目标非线性规划问题的求解2.1 收敛性分析针对问题P1, 构造如下动态系统的状态方程:

(15)

(16)

(17)

(18)

当动态系统处于平衡点时, 系统内各变量(x, s, t, q)处于稳定状态, 各变量的时间导数为0, 可得如下关系式:

(19)

(20)

(21)

(22)

其中, s和t为约束条件相对应的非负松弛变量, 则对于任意的x, 必有一组s和t使得式(20)和(21)成立.设决策变量x和松驰变量q的初始值分别为x₀和q₀, 则初始条件下问题P1的目标函数值为f(x₀), 此时f(x₀)大于最优解的目标函数值f(x^*).在动态系统存在稳定状态的情况下, 系统中各变量随时间达到稳定状态.根据式(19), f (x) 和q的取值将逐渐接近直至相等.
根据松弛变量q和目标函数值f (x) 在系统收敛过程中的变化情形来对李雅普诺夫方法求解的质量进行判断.根据目标函数的松弛变量初值q₀与f(x^*) 的大小关系、决策变量x和松弛变量q的变化速度的不同, 将动态系统的收敛情形主要分为以下三种:
情形1 q₀≥f(x^*).
系统的收敛过程如图 1所示.此时存在f(x₀) ≥q₀和f(x₀) < q₀两种情况, 因为f(x₀)≥ f(x^*)成立, 所以, 无论x和q在系统收敛过程中的变化速度如何, f(x)和q在逐渐汇合的过程中, 二者的值始终会大于f(x^*).稳定状态时得到的目标函数值f(x_end)必然大于f(x^*).因此, 此时得到的解x_end不保证是问题P1的最优解.
图 1(Fig. 1)

图 1 情形1下系统的收敛过程Fig.1 Convergence process of the system in situation 1 (a)—f(x0) ≥q0情形；(b)—f(x0) <q0情形.

情形2 q₀ < f(x^*) 且q的变化速度较快.
因为松弛变量q的变化速度较快, 所以将在系统收敛过程中率先到达f(x^*) 处, 系统的收敛过程如图 2所示.由于松弛变量q的变化速度较快, 为满足式(19)成立, q会越过f(x^*) 的位置, 直至与f(x) 汇合.最终得到的解x_end对应的目标函数值f(x_end)将大于f(x^*), 所以此时求得的解也不保证是问题P1的最优解.
图 2(Fig. 2)

图 2 情形2下系统的收敛过程Fig.2 Convergence process of the system in situation 2

情形3 q₀ < f(x^*) 且q的变化速度较慢.
如图 3所示, 因为松弛变量的变化速度较慢, f(x) 在收敛过程中将率先到达f(x^*) 处.与情形2相似, 为了满足式(19)成立, f(x) 会越过f(x^*) 处向q靠近.但是, 由于f(x) 受到约束条件的限制, f( x) 将难以继续向q靠近.在与松弛变量q汇合前, f(x) 将在靠近f(x^*) 处震荡变化, 且幅度越来越小, 直到与q汇合, 此时可认为f(x) 与f(x^*) 相等.因此, 此情形下松弛变量q与f(x) 汇合得到的解x_end是问题P1一定误差范围内的最优解.
图 3(Fig. 3)

图 3 情形3下系统的收敛过程Fig.3 Convergence process of the system in situation 3

综上所述, 情形3可保证动态系统收敛到最优解.此时应满足如下条件: 松弛变量q的初始值q₀ < f(x^*), 且在收敛过程中决策变量x的目标函数值f (x)先于q到达f (x^*) 处.
2.2 参数设置基于上述分析, 在应用李雅普诺夫方法时, 给出如下参数设置的建议:
1) k_q? k_x, 其中“?”表示“远小于”.此时决策变量 x变化速度将远快于松弛变量q, 从而可以保证f(x) 更快到达最优解f(x^*) 处, 可避免情形2的出现.
2) q₀?f(x^*).此时, f(x^*)- q₀?f(x₀)- f(x^*) 成立.由于松弛变量的初始值q₀到f(x^*) 的距离远大于f(x₀) 到f(x^*) 的距离, 即使k_x和k_q不严格满足1)的建议, 仍可保证f(x)率先到达f(x^*) 处, 可避免情形1的出现, 保证求解质量.
2.3 算例分析通过实验对上述结论进行说明性验证, 求解工具为MATLAB软件的ODE求解器.待求解的约束规划问题为(例1)：

(23)

(24)

已知最优解为x₁^*=1.105 6, x₂^*=0.552 9, 目标函数值f(x^*) =1.527 9.
2.3.1 三种情形的验证分3种情形设置参数如下: 情形1, q₀=3, k_x=(1，1)，k_q=1;情形2, q₀=0, k_x=(0.1，0.1), k_q=1;情形3, q₀=0, k_x=(10，10), k_q=1. 三种收敛情形下决策变量x对应的目标函数值f (x) 和松弛变量值q随仿真时间的变化情况如表 1所示.
表 1(Table 1)

表 1 3种情形下目标函数值和松弛变量随时间的变化情况Table 1 The change of objective function value and slack variable with time in three situations 收敛情形 仿真时间/s 目标函数f(x) 松弛变量q 情形1 10-6 2.000 1 3.000 0 情形1 10-5 2.000 8 3.000 0 情形1 10-4 2.008 0 3.000 0 情形1 10-3 2.077 6 2.999 0 情形1 10-2 2.570 3 2.993 3 情形1 10-1 2.981 6 2.987 3 情形1 10-0 2.987 1 2.987 1 情形1 101 2.987 1 2.987 1 情形1 102 2.987 1 2.987 1 情形2 10-6 2.000 0 0.000 2 情形2 10-5 1.999 8 0.002 0 情形2 10-4 1.998 4 0.020 0 情形2 10-3 1.984 9 0.189 6 情形2 10-2 1.905 6 1.224 0 情形2 10-1 1.878 8 1.876 7 情形2 100 1.928 1 1.928 0 情形2 101 1.929 7 1.929 7 情形2 102 1.929 7 1.929 7 情形3 10-6 1.999 8 0.000 0 情形3 10-5 1.998 4 0.000 0 情形3 10-4 1.984 1 0.000 2 情形3 10-3 1.853 4 0.001 9 情形3 10-2 1.208 3 0.015 0 情形3 10-1 0.922 7 0.095 1 情形3 10-0 1.103 9 0.650 6 情形3 101 1.516 3 1.508 7 情形3 102 1.527 9 1.527 9

表 1 3种情形下目标函数值和松弛变量随时间的变化情况 Table 1 The change of objective function value and slack variable with time in three situations

为方便分析, 将3种情形的收敛结果展示在图 4~图 6中.
图 4(Fig. 4)

图 4 例1情形1下求解结果Fig.4 Results of example 1 in situation 1

图 5(Fig. 5)

图 5 例1情形2下求解结果Fig.5 Results of example 1 in situation 2

图 6(Fig. 6)

图 6 例1情形3下求解结果Fig.6 Results of example 1 in situation 3

由图 4和5可知, 虽然目标函数值f (x) 和松弛变量值q相互接近直至相等, 但二者的重合点在f(x^*) 之上, 即f (x_end)>f (x^*), 分别对应于情形1和情形2.由图 6可知, 在动态系统收敛的初始阶段, f(x) 变化速度较快, f(x) 迅速向q靠近并率先抵达f(x^*) 处, 但在越过f(x^*) 处后, 其变化速度明显下降, 直至反向变化, 最终稳定在f(x^*) 处.松弛变量q则是一直向f(x) 方向变化最终与之重合.这对应于情形3, 可以得到问题的最优解.
2.3.2 参数选择与求解质量为验证上述参数设定的建议, 改变k_x, k_q和q₀的取值, 重新求解, 所得结果如表 2所示.可见, 在保证q₀ < f(x^*) 的前提下, 松弛变量初值q₀越接近最优目标函数值f(x^*), 其对k_x, k_q间的大小差距要求越大；反之, 则对k_x, k_q间的差距要求越低.选取q₀=1.5, 1-1，1-5，1-9，1-13实验为一组, 此时q₀已接近f(x^*), 只有k_x/k_q=100才能保证理想的求解精度, 而其余实验无法求得最优解, 且随着k_x/k_q的减小, 求解质量也有明显下降.
表 2(Table 2)

表 2 例1在不同参数下的求解结果Table 2 Results of example 1 under different parameters 序号 kx/kq q0 求解结果xend 目标值f(xend) 误差/% 1-1 100 1.5 (1.105 6, 0.552 9) 1.527 9 0.00 1-2 100 0.0 (1.105 6, 0.552 9) 1.527 9 0.00 1-3 100 -1.5 (1.105 6, 0.552 9) 1.527 9 0.00 1-4 100 -3.0 (1.105 6, 0.552 9) 1.527 9 0.00 1-5 10 1.5 (1.073 5, 0.623 6) 1.541 4 0.88 1-6 10 0.0 (1.105 6, 0.552 9) 1.527 9 0.00 1-7 10 -1.5 (1.105 6, 0.552 9) 1.527 9 0.00 1-8 10 -3.0 (1.105 6, 0.552 9) 1.527 9 0.00 1-9 1 1.5 (1.024 9, 0.778 0) 1.655 8 8.37 1-10 1 0.0 (1.105 0, 0.553 8) 1.527 9 0.00 1-11 1 -1.5 (1.105 6, 0.552 9) 1.527 9 0.00 1-12 1 -3.0 (1.105 6, 0.552 9) 1.527 9 0.00 1-13 0.1 1.5 (1.002 3, 0.932 1) 1.873 5 22.62 1-14 0.1 0.0 (1.025 3, 0.776 5) 1.654 2 8.27 1-15 0.1 -1.5 (1.051 0, 0.684 8) 1.573 5 2.98 1-16 0.1 -3.0 (1.070 0, 0.633 0) 1.545 1 1.13

表 2 例1在不同参数下的求解结果 Table 2 Results of example 1 under different parameters

然后, 选取k_x/k_q=1, 1-9，1-10，1-11，1-12实验为一组, 可以发现实验1-9未求得最优解, 实验1-10得到了一定误差范围内的最优解, 实验1-11，1-12则求得了最优解.这说明随着q₀与f(x^*) 大小差距的增大, 算法的求解质量在逐渐提升.这与上文的分析一致, 从而验证了前文论述的正确性.
综上, 依据2.2节所述建议对参数进行设定可以保证求解结果的最优性, 但是不建议为了保证求解质量盲目增大q₀与f( x^*) 的差距及k_x/k_q的大小.在一定范围内增加k_x/k_q值可以提高求解精度, 但超过一定范围后, 将无法再明显提高精度, 反而会使系统收敛困难, 增加求解时间; 同理, 虽然增大q₀与f(x^*) 的差距可以保证求解质量, 但在超过一定范围后, 将无法再对求解质量产生明显影响^[14-15].根据实验1-5，1-6，1-7，1-8求解结果可验证以上结论.因此, 在对李雅普诺夫方法的初始参数设定时, 需要考虑问题的实际情况和精度要求.
下面考虑约束条件相对应的非负松弛变量s和t对变量x的影响, 待求解的约束规划问题为(例2):

(25)

(26)

已知最优解x^*=3, 目标函数值f(x^*)=-9.改变非负松弛变量s和t, 构造方程组并求解, 所得结果如表 3所示.
表 3(Table 3)

表 3 例2在改变参数s和t时的求解结果Table 3 Results of example 2 when the parameters s and t are changed 序号 s t 求解结果x 2-1 1 1 3 2-2 10 1 3 2-3 100 1 3 2-4 200 1 3 2-5 1 10 3 2-6 1 100 3 2-7 1 200 3

表 3 例2在改变参数s和t时的求解结果 Table 3 Results of example 2 when the parameters s and t are changed

从表 3可以看出, 非负松弛变量s和t的数值选取不会影响最优解的求解.
3 多目标非线性规划问题的求解对于多目标非线性规划问题的求解, 元启发式算法(如遗传算法^[16]、进化规划算法^[17]、粒子群优化法^[18]等)虽可以取得较好的效果, 但是多数情况下收敛性都悬而未决^[14].本节分析应用李雅普诺夫方法解决多目标非线性规划问题的收敛性及参数设置.
3.1 收敛性分析给出多目标优化问题的一般形式(P2):

(27)

(28)

(29)

其中：f(x) 为p维目标函数；g(x) 为n维等式约束；h(x) 为m维不等式约束, h_min和 h_max分别为不等式约束的下界和上界.
针对多目标规划问题P2, 构造动态系统, 得到系统的状态方程为

(30)

(31)

(32)

(33)

当该系统收敛到稳定状态时, 各变量的状态方程满足, 由此可得如下关系式:

(34)

(35)

(36)

(37)

(38)

此时求得的解是多目标优化问题P2的可行解.
设决策变量x和松弛变量q的初始值分别为x₀和q₀, 则初始条件下问题P2目标函数值为f(x₀), 此时f(x₀) 位于原问题可行解的目标函数f(x) 的值域内.与单目标规划问题P1相似, 通过分析动态系统在不同初始条件下的收敛情形, 得到该方法在解决多目标规划问题时不同的求解结果.根据松弛变量初值q₀与目标函数值域的位置关系, 以及决策变量x和松弛变量q的变化速度的关系, 将动态系统的收敛情形主要分为以下三种:
情形1 q₀在目标空间内.
q₀位于问题P2的目标空间之内, 收敛过程如图 7所示.因为f(x₀) 处于目标函数值的值域之内, 所以无论x和q在系统收敛过程中的变化速度如何, 在f(x) 和q逐渐靠近的过程中, 二者数值会一直处于目标函数值域之内.因此, f(x) 和q最终汇合得到的解对应的目标函数f(x_end)不保证位于目标空间的边界上, 此时求得的解x_end不保证是问题P2的Pareto解.
图 7(Fig. 7)

图 7 多目标问题情形1下的收敛情形Fig.7 Convergence process of multi-objective problems in situation 1

情形2 q₀位于目标空间之外且松弛变量的变化速度较快.
q在系统收敛过程中率先到达问题P2的目标空间边界处, 收敛情形如图 8所示.由于松弛变量q会越过目标空间边界继续向f(x) 的方向变化, 直至二者汇合.此时最终得到的目标函数值f(x_end)将会处于问题P2目标空间内而不是边界处, 所以此时求得的解不保证是问题P2的Pareto解.
图 8(Fig. 8)

图 8 多目标问题情形2下的收敛情形Fig.8 Convergence process of multi-objective problems in situation 2

情形3 q₀位于目标空间之外且松弛变量的变化速度较慢.
f(x) 在收敛过程中率先到达目标空间的边界处, 收敛过程如图 9所示. f(x) 由初始位置出发向q₀靠近的过程中由于受约束条件的限制无法越过目标空间, 所以这一阶段f(x) 只能停留在该空间的边界上, 且距离q₀最近; 同时, q与f(x^*) 逐渐汇合的过程中, q值的变化可能会影响到f(x), 表现为f(x) 在目标空间的边界上滑动.因此, 此情形下只要通过调整q₀与目标函数值域的位置关系, 令其位于值域的左下方, 此时f(x) 就会在收敛过程中向目标函数值域的左下边界靠近并停留.此时松弛变量 q与f(x) 汇合得到的解x_end是问题P2一定误差范围内的Pareto最优解.
图 9(Fig. 9)

图 9 多目标问题情形3下的收敛情形Fig.9 Convergence process of multi-objective problems in situation 3

3.2 参数设置由上述分析可知, 为了保证应用李雅普诺夫方法求解多目标规划问题P2的结果为Pareto最优解, 需要保证两个条件: 1) q₀位于问题P2的目标空间之外, 保证q₀小于目标函数下界; 2)在动态系统收敛的过程中相对于松弛变量q, f(x) 先抵达问题的目标函数值域边界处.当条件1)中目标空间无法确定时, 可以令q₀的值足够小, 无限接近无穷小, 此时f(x) 若能以足够快的速度向目标值f(x^*) 移动, 算法仍可以使用.
不同的参数选择也会影响动态系统最终的平衡位置, 从而得到不同的Pareto解.下面以2个目标函数f_i(i=1, 2)的规划问题为例, 讨论松弛变量q_0i、松弛变量对应的增益因子k_qi对动态系统平衡点的影响.
1) 松弛变量的增益因子k_qi.假设q₁的收敛速率比q₂慢, 则在收敛的过程中, q₂就会在f₂的方向上移动更大的距离.这也会导致在动态系统收敛的前一阶段, f(x) 在问题目标空间边界上的停留位置容易向f₁小的方向靠拢.因此, 与其他松弛变量的增益因子相比较, k_qi值越大, 最终求得的解对应的目标值中f_i^*则会越大.
2) 松弛变量的初值q_0i.假设q₀₂到问题P2的目标空间的距离比q₀₁更远, 那么f(x) 与q收敛过程中, f(x) 会在距离相对较大的f₁方向上移动更大的距离, f(x) 最终在目标空间边界停留的位置也会向f₁较小的方向靠近.因此, 松弛变量的初值q_0i越小, 最终求得的解对应的目标值f_i^*则会越小; 反之, q_0i越大, 最终求得的解对应的目标值f_i^*则会越大.
3.3 算例分析

(39)

(40)

(41)

(42)

(43)

(44)

式(39)~(44)为多目标规划问题的一个例子(例3).不同参数下应用李雅普诺夫方法求解的结果如表 4所示.
表 4(Table 4)

表 4 例3在不同参数下的求解结果Table 4 Results of example 3 under different parameters 序号 增益因子(kq1, kq2) 松弛变量(q01, q02) x1 x2 f1(xend) f2(xend) 3-1 (10, 1) (-500, -500) -14.981 9 0.736 2 290.455 3 -134.768 0 3-2 (8, 1) (-500, -500) -14.982 4 0.725 4 290.478 9 -134.767 0 3-3 (6, 1) (-500, -500) -14.983 5 0.704 3 290.525 2 -134.764 0 3-4 (4, 1) (-500, -500) -14.131 0 0.999 9 262.209 1 -127.179 0 3-5 (2, 1) (-500, -500) -8.279 5 1.000 0 107.668 7 -74.515 7 3-6 (1, 1) (-300, -300) -2.119 4 2.626 9 21.616 0 -16.427 8 3-7 (1, 1) (-300, -600) -5.991 6 1.633 4 52.799 5 -45.490 9 3-8 (1, 1) (-300, -900) -7.532 4 1.000 0 92.867 4 -67.792 0 3-9 (1, 1) (-300, -1 200) -9.528 9 1.000 0 134.916 7 -85.760 5 3-10 (1, 1) (-300, -1 500) -11.211 6 1.000 0 176.547 6 -100.905 0

表 4 例3在不同参数下的求解结果 Table 4 Results of example 3 under different parameters

以序号3-1到3-5五个实验为一组, 观察k_q不同取值下的求解结果, 可以得出：随着k_q1减小, 目标函数值f₁(x_end)逐渐减小, 但f₂(x_end)逐渐增大; 以序号3-6，3-8，3-9和3-10为一组, 观察改变松弛变量q₀的求解结果, 可以得出随着q₀₂减小, 目标函数值f₂(x_end)也在逐渐减小, 但f₁(x_end)逐渐增大, 从而验证了3.2节所得的结论.
进而, 改变f₁, f₂的权重, 采用加权法对例3问题进行求解, 结果如表 5所示.两种方法求解结果的对比如图 10所示.可见: 利用加权法, 在不同权值下求得的解分布比较集中, 可以看作是在一定误差范围内应用李雅普诺夫方法求解多目标问题的所得解的特例.相比之下, 应用李雅普诺夫方法求解能够给出满足约束条件的非劣解集, 更方便决策者根据实际情形进行决策.因此, 针对解决多目标规划问题, 应用李雅普诺夫方法比加权法更加有效.
表 5(Table 5)

表 5 不同权重比例下例3的求解结果Table 5 Results of example 3 with different weights ratio 序号 权重比例 x1 x2 目标值f1(xend) 目标值f2(xend) 4-1 (0.1, 0.9) 0 3.333 3 11.444 3 5.444 3 4-2 (0.2, 0.8) 0 3.333 3 11.444 3 5.444 3 4-3 (0.3, 0.7) 0 3.333 3 11.444 3 5.444 3 4-4 (0.4, 0.6) 0 3.333 3 11.444 3 5.444 3 4-5 (0.5, 0.5) 0 3.333 3 11.444 3 5.444 3 4-6 (0.6, 0.4) 0 3.333 3 11.444 3 5.444 3 4-7 (0.7, 0.3) 0 3.333 3 11.444 3 5.444 3 4-8 (0.8, 0.2) 0 3.333 3 11.444 3 5.444 3 4-9 (0.9, 0.1) 0.56 3.333 3 10.417 1 11.453 6

表 5 不同权重比例下例3的求解结果 Table 5 Results of example 3 with different weights ratio

图 10(Fig. 10)

图 10 李雅普诺夫方法与加权法求的对比Fig.10 Comparison between Lyapunov method and weighted method

4 结语本文分别针对李雅普诺夫方法求解单目标和多目标约束非线性规划问题的求解过程进行分析, 给出了算法的收敛性, 并提出算法中关键参数的取值建议.对于最小化单目标非线性规划问题, 当松弛变量初始值q₀小于问题的最优值, 且在收敛过程中决策变量x比松弛变量q变化速度更快时, 可保证系统收敛到最优解, 松弛变量s和t的取值不影响最优解的求解; 对于最小化多目标约束非线性规划问题, 松弛变量的增益因子k_qi和松弛变量的初值q_0i越小, 对应的目标值f_i则会越小.最后给出了数值算例.
李雅普诺夫方法是求解约束非线性规划问题的一种新颖的方法.未来可以对该算法的并行式实施、分布式实施进行研究, 与现有算法进行更大范围的对比, 并尝试应用于求解智能制造、边缘计算等领域更多的实际问题, 这些也是作者未来的研究方向.
参考文献

[1]	Fliege J, Drummond L M G, Svaiter B F. Newton's method for multiobjective optimization[J]. SIAM Journal on Optimization, 2009, 20(2): 602-626. DOI:10.1137/08071692X
[2]	Rao R V, Savsani V J, Vakharia D P. Teaching-learning-based optimization: an optimization method for continuous non-linear large scale problems[J]. Information Sciences, 2012, 183(1): 1-15. DOI:10.1016/j.ins.2011.08.006
[3]	Hu J, Jiang B, Lin L, et al. Structured quasi-Newton methods for optimization with orthogonality constraints[J]. SIAM Journal on Scientific Computing, 2019, 41(4): A2239-A2269. DOI:10.1137/18M121112X
[4]	Liu X, Yuan Y. A sequential quadratic programming method without a penalty function or a filter for nonlinear equality constrained optimization[J]. SIAM Journal on Optimization, 2011, 21(2): 545-571. DOI:10.1137/080739884
[5]	Hong M, Luo Z Q, Razaviyayn M. Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems[J]. SIAM Journal on Optimization, 2016, 26(1): 337-364. DOI:10.1137/140990309
[6]	Li G, Pong T K. Global convergence of splitting methods for nonconvex composite optimization[J]. SIAM Journal on Optimization, 2015, 25(4): 2434-2460. DOI:10.1137/140998135
[7]	Han C, Wang L, Zhang Z L, et al. A multi-objective genetic algorithm based on fitting and interpolation[J]. IEEE Access, 2018, 6(1): 22920-22929.
[8]	Chen Y H, Zhang B. Multi-objective optimization of current limiting scheme based on NSGA-Ⅱ[C]// Joint Conference of 6th International Congress on Electrical and Engineering(ICECE). Shanghai, 2019: 205-214.
[9]	Gadhvi B, Savsani V, Patel V. Multi-objective optimization of vehicle passive suspension system using NSGA-Ⅱ, SPEA2 and PESA-Ⅱ[C]//3rd International Conference on Innovations in Automation and Mechatronics Engineering. Vallabh Vidhyanagar, 2016: 361-368.
[10]	Torelli F, Vaccaro A, Xie N. A novel optimal power flow formulation based on the Lyapunov theory[J]. IEEE Transactions on Power Systems, 2013, 28(4): 4405-4415. DOI:10.1109/TPWRS.2013.2266126
[11]	Torelli F, Montegiglio P, De Bonis A, et al. A new approach for solving DAE systems applied to distribution networks[C]// 49th International Universities Power Engineering Conference(UPEC). Cluj-Napoca, 2014: 131-139.
[12]	Torelli F, Vaccaro A. A generalized computing paradigm based on artificial dynamic models for mathematical programming[J]. Soft Computing, 2014, 18(8): 1561-1573. DOI:10.1007/s00500-013-1162-z
[13]	Torelli F, Vaccaro A. A second order dynamic power flow model[J]. Electric Power Systems Research, 2015, 126(1): 12-20.
[14]	Liang X, Mosalam K M. Lyapunov-based nonlinear solution algorithm for structural analysis[J]. Journal of Engineering Mechanics, 2018, 144(9): 12-21.
[15]	Chellaboina V, Leonessa A, Haddad W M. Generalized Lyapunov and invariant set theorems for nonlinear dynamical systems[J]. Systems & Control Letters, 1999, 38(45): 289-295.
[16]	Xie N, Bompard E, Napoli R, et al. Widely convergent method for finding solutions of simultaneous nonlinear equations[J]. Electric Power Systems Research, 2012, 83(1): 9-18. DOI:10.1016/j.epsr.2011.09.002
[17]	Xie N, Torelli F, Bompard E, et al. Dynamic computing paradigm for comprehensive power flow analysis[J]. IET Generation Transmission & Distribution, 2013, 7(8): 832-842.
[18]	Fonseca C M, Fleming P J. Multiobjective optimization and multiple constraint handling with evolutionary algorithms Part I: unified formulation[J]. IEEE Transactions on Systems Man and Cybernetics Part A-Systems and Humans, 1998, 28(1): 26-37. DOI:10.1109/3468.650319