第二章 鲁棒控制理论概述

2．1鲁棒控制理论概述 2．1．1 系统不确定性和鲁棒性

控制科学所要解决的主要问题之一是针对被控对象，设计合适的控制器，使闭环系统稳定或达到一定的性能指标要求。它经历了经典控制理论和现代控制理论两个发展阶段。无论是经典控制理论还是现代控制理论，它们的一个明显的特点是建立在精确的数学模型基础之上。但是，在实际应用中存在着许多不确定性，具体体现在：

(1)参数的测量误差。由于测量技术的，许多参数的测量值可能有相当大的误差。尤其是某些涉及热力学、流体力学和空气动力学，以及化学反应过程的参数，往往很不容易测准，或者需要付出昂贵的代价才能测准；

(2)环境和运行条件的变化。这往往是不确定性产生的最重要的原因。例如，内部元器件的老化；电气设备的电阻因温升而改变；炼钢炉因炉壁渐渐被钢水腐蚀变薄而导致导热系统的变化；飞机和导弹在高空或低空以高速或低速飞行时其空气动力学参数的变化非常剧烈，甚至由于燃料消耗造成导弹质量的变化和质心的位移，这些都会造成其参数较大的变化； (3)人为的简化。为了便于研究和设计，人们往往有意略去系统中一些次要因素，用低阶的线性定常集中参数模型来代替实际的高阶、非线性甚至是时变和分布参数的系统，这样势必要引入系统模型的不确定性。因此，在控制系统的设计过程中不可避免的问题是：如何设计控制器，使得当一定范围的参数不确定性及一定限度的未建模动态存在时，闭环系统仍能保持稳定并保证一定的动态性能，这样的系统被称为具有鲁棒性。

2．1．2鲁棒控制理论的发展概况

鲁棒控制理论正是研究系统存在不确定性时如何设计控制器使闭环系统稳定且满足一定的动态性能。自从1972年鲁棒控制(Robust Contr01)这一术语首次在期刊论文中出现以来，已有大量的书籍详细的阐述了鲁棒控制理论的产生、发展及研究现状。鲁棒控制的早期研究常只限于微摄动的不确定性，都是一种无穷小分析的思想。1972年鲁棒控制(Robust Control)这一术语首次在期刊论文中出现。经过三十多年的研究，鲁棒控制理论已比较成熟，在时域和频域都取得了令人瞩目的成就，其代表性的研究方法有多项式代数方法以Kharitonov定理为代表的多项式代数方法，为参数不确定系统的鲁棒控制研究提供了强有力的理论方法，但由于本身理论的局限性，此方法基本上只能局限于多项式空间和对系统鲁棒稳定性的分析，对参数不确定系统的鲁棒镇定问题，一直没有什么满意的结果。如何将现有方法应用到控制工程实践，仍有许多问题需要解决。H控制理论的提出具有很强的工程应用背景。控制理论的基干扰信号属于某一有限能量信号集情况下，用其相应的灵敏度函数标，从而将干扰问题化为求解使闭环系统稳定，并使相应的如范数馈控制问题。比设计方法虽然将鲁棒性直接反映在系统的设计指标映在相应的加权函数上，但它“最坏情况”下的控制却导致了

1

不必要的保守性。另外，由于巩设计方法是以非结构化不确定性和小增益定理为设计框架的由于巩设计方法是以非结构化不确定性和小增益定理为设计框架的，稳定性，对鲁棒性能的分析就显得为力。因此，鲁棒多变量反馈系统设计方法一直存在的困难，是不能够在统一的框架下同时处理性能指标与鲁棒稳定性的折中问题。

μ方法很好地补充了比控制的不足，可以把结构不确定系统的鲁棒性能结合起来考虑，并且克服了设计上的保守性，从而设计出性能更优，鲁棒性更好的控制系统。值得注意的是，μ理论中的μ分析已基本完善，但μ综创还没有很好地解决。目前常用的是Doyle提出的D一Ⅳ迭代算法，由于Ⅳ和D的优化并不具有组合凸性，所以不能保证迭代算法收敛到全局最优，因而求得的弘值具有一定的保守性，在一定程度上了“理论的具体应用。鲁棒控制理论的时域法是鲁棒控制理论中最活跃的分支，它考虑实际系统与数学模型之间存在偏差时，如何保证系统的稳定性和其他性能。自从Monopolill首次采用Lyapunov稳定性理论研究不确定系统的鲁棒镇定问题以来，对于时变和非线性摄动不确定系统，基于Lyapunov稳定性理论的鲁棒镇定综合方法引起了众多学者的关注。在这一框架内主要有两种研究方法，即Riccati方程处理方法和线性矩阵不等式方法。Riccati方法是早期的一种研究方法，其基本思想是将不确定系统的分析和综合问题转化为一个Riccati方程(或不等式)的可解性问题，进而通过求解Riccati方程来对系统的鲁棒稳定性和鲁棒性能进行分析，或给出鲁棒控制器。Riccati方程处理方法在80年代和90年代初期被广大学者采用，但随着研究问题的日益复杂，越来越多的学者认识到Riccati方法的局限性：

(1)Riccati方程的求解存在一定的问题。虽然目前有很多求解Riccati方程的方法，但大多为迭代方法，其收敛性不能得到保证；

(2)众所周知，应用Riccati方法进行不确定系统的分析和综合时，往往需要设计者预先确定一些待定参数，这些参数的选择不仅直接影响到结论的好坏，而且还会影响到问题的可解性。在现有的Riccati方程处理方法中，还缺乏寻找这些参数最佳值的方法，多数情况下尚需要 人为的确定这些参数，无疑给分析和综合结果引入很大的保守性。随着求解凸优化问题的内点法的提出，到20世纪90年代初，线性矩阵不等式方法逐渐受到控制界的普遍关注。通过线性矩阵不等式技术，系统和控制中的很多问题可以转化为一个线性矩阵不等式(组)的可行性问题，或者转化为一个受线性矩阵不等式(组)约束的凸优化问题。内点法的提出使鲁棒分析和综合中的一些原来无法解决的复杂问题在转化为线性矩阵不等式问题后得以有效的解决。尤其是MathWork公司在其商业软件Matlab中推出了求解线性矩阵不等式问题的LMI工具箱后，使得人们在求解线性矩阵不等式问题时更方便、更有效，这进一步推动了线性矩阵不等式方法在系统和控制领域中的应用。线性矩阵不等式处理方法克服了Riccati方程方法中的许多不足。采用线性矩阵不等式方法处理不确定系统的鲁棒分析和综合问题时，所需要预先选择的参数要明显少于Riccati方法。线性矩阵不等式方法给出了问题解的一个凸约束条件，可以采用求解凸优化问题的有效方法，得到一组满足要求的可行解，而不是唯一 解，因而可以对这一组解做进一步优化，这也使得线性矩阵不等式方法不仅为广大科研工作者所采用，而且正逐渐为工程应用人员所接纳。

2．2 时滞系统鲁棒控制概述

动力系统总是存在滞后现象。从工程技术、物理、力学、控制论、化学反应、生物医学等中提出的数学模型带有明显的滞后量，特别是在自动控制的装置中，任何一个含有反馈的系统，从输入信号到收到反馈信号，其间必然有一个时间差。因此时滞是普遍存在的。例如在化工、液压、轧钢、核反应堆、轮船定向仪、无损传输系统等系统中都具有时滞，而且时滞是引起系统不稳定以及导致系统性能恶化的一个重要因素。因此，对时滞系统的研究具有重要的理

2

论意义与应用价值。 考虑线性时滞系统

x(t)Ax(t)Adx(xd)x(t)(t)dx0 （2.2.1）

nnn其中x(t)R为状态变量，d≥0为时滞，(t)为初始条件，矩阵A,AdR为知的定常

矩阵。早期对时滞系统的研究通常采用频域的方法。一般时滞系统有无穷多个极点，很难用传统的方法将这些极点配置到左半复平面的指定位置。Smith预估是克服难的一种方法。状态预估器和过程模型控制采用了Smith预估器的思想，解决了指定干扰抑制问题。但是，以上方法仅限于标称对象的控制问题。当系统存在不确定性时，这些方法就为力。近年来，许多学者做了改善研究。然而，对于时变时滞还是未能得到有效地解决。从时域的角度，Lyapunov泛函方法是处理时滞系统的一般方法，其主要思想是通过构造一个合适的Lyapunov-Krasovskii泛函或Lyapunov函数，获得系统（2.2.1）稳定的充分条件。它克服了频域不能处理时变和参数摄动的不足，而且随着线性矩阵不等式技术的成熟，使其计算简单，因而在时滞系统的分析和设计中得到了广泛地应用。 在20世纪90年代前，提出的关于时滞系统的结论基本上都是时滞无关的，也就是说在分析 和设计系统时，不考虑系统实际时滞的大小，因而所得结论对任意大小的时滞都成立，这对于无法精确得到系统滞后信息的一类时滞系统无疑是有效的，但当时滞很小时，这种时滞无关(delay-independent)条件是相当保守的。相对应地，当考虑时滞对系统性能的影响时得出的条件就称为时滞相关(delay-dependent)条件。这类条件须首先假设当d=0时系统(1．1)是稳定的，这样由于系统的解对d的连续依赖，则一定存在一个时滞上界dm∞，使得系统（2.1）对dx0都是稳定的。因此，最大容许的时滞上界就成为衡量时滞相关条件保守性的主要指标。近年来，时滞相关稳定性分析与控制综合，以及如何降低所得条件的保守性，已经成为控制理论界研究的热点问题。目前，国际上主要针对两类时滞研究其时滞相关问题，第一类是定常时滞、第二类是时变时滞。而时滞时变情形又分为两种情形(Case I)和(CaseⅡ)，其中Case I是指时变时滞连续但不可微，Case U是指时变时滞连续且可微。在研究方法上，不管是定常时滞还是时变时滞，主要采用的是时域研究方法，可分为四类：离散Lyapunov泛函方法、模型变换方法、自由权矩阵方法、积分不等式(有限和不等式)方法。离散Lyapunov泛函方法的基本思想是对Lyapunov泛函进行离散化，获得LMI表示的系统稳定性结果。该方法的优点是：只要步长足够小，对于保证系统稳定的时滞界限的估计就非常接近于实际值，但其算法复杂，而且该方法只适用于具有定常时滞的系统，对于具有时变时滞的系统就为力；另外，该方法只适用于稳定性分析，很难推广到综合问题。因此，这类方法自从1997年Gu提出后，只有少数学者进行研究，没有得到广泛的推广和应用。模型变换方法主要是将一个具有离散时滞的系统通过Leibniz．Newton公式, 将线性时滞系统（2.1）转化为一个具有分布时滞的新系统，再选取适当的Lyapunov泛函，从而得到时滞相关条件。模型变换的目的是让系统方程中出现积分项，这样对y函数沿系统求导就导致交叉项与二次型积分项的同时出现，然而对交差项的界定可以抵消y泛函导数中的二次型积分项，从而可获得时滞相关条件。这3种模型变换方法简单，对于稳定性和性能分析，基本上都能用LMI求解。而且能够推广考虑各种综合问题来求解控制器。特别是基于Park不等式和Moon不等式来界定交叉项的模型变换，已经有一系列的结果，如：有界实引理，H控制，、鲁棒稳定性等。但是这些守性，一方面是模型变换带来了保守性；另一方面对交叉项的界定也不可避免地产生保守性。

3

2.3 LMI设计实例：风力发电机整流器设计

在风力发电控制系统中，变流器是接在发电机和电网之间的。而网侧变流器（三相PWM整流器）在工作时，能够在稳定直流侧电压的同时，实现其交流侧在受控功率因数（如单位功率因数）条件下的正弦波电流控制。另外一方面，常规的三相电压型PWM整流器（VSR）控制系统通常应用双闭环的控制策略，即电压外环控制和电流内环控制。

2.3.1、电流环的名义系统模型

固定开关频率PWM电流控制算法简单、滤波电感的设计比较容易，且实现起来比较为方便。所以本设计中采用直接电流控制下的固定开关频率PWM电流控制。一般我们把固定开关频率PWM电流控制方案又分为静止abc坐标系下的电流控制方案和旋转dq坐标系下的电流控制方案。三相电压型PWM整流器在d,q同步旋转坐标系下的dq模型也可以描述为

edLpRLidvdeiv (2.3.1)

LLpRqqq3(vdidvqiq)vdcidc (2.3.2) 2式（1）中ed、eq为电网电动势矢量Edq在d、q轴的分量；vd、vq是三相VSR交流侧电压矢量Vdq的在d、q轴的分量；id、iq是三相VSR交流侧电流矢量Idq在d、q轴的分量；p为微分算子。

假设d,q同步旋转坐标系下的q轴与电网电动势矢量Edq是重合的，那么电网电动是矢量在d轴的分量ed为0。从式（1）中可以看出，由于三相电压型PWM整流器在d,q轴的变量是相互耦合的，所以在控制器设计时就形成了困难。那么，我们选用前馈解耦的控制方案，电流内环采用PI调节器，因此vd和vq的控制方程为

vq(KiPKiI*)(iqiq)Lideq (2.3.3) sK*vd(KiPiI)(idid)Liqed (2.3.4)

s*上式中，KiI和KiP是电流环的积分调节增益和比例调节增益；而id、iq为id、iq的电流指令。我们把式（3）和式（4）代入式（1）后，可得

* 4

KiI[R(K)]/L0*iPid1idKiIidsp)*（2.3.5) (KiPiiKLsqiq0[R(KiPiI)]/Lqs从式（5）我们可以很明显的看出，式（3）和式（4）使三相电压型PWM整流器的电流环实现了解耦。因为两电流环的对称性，我们仅以iq的控制为例来设计电流调节器。这里考虑PWM控制的小惯性特性和电流环的信号采样的延迟。T是电流环的电流采样周期（即PWM开关周期），桥路PWM的等效增益是K。为了方便计算，eq的扰动这里暂不考虑，忽略负载电流iL的扰动，那么，三相电压型PWM整流器的电流环控制结构图见图1

Iq_in Iq_out 1 Ts1K 0.5Ts11/R L/Rs1图1电流环控制结构框图

从结构图1可以得到电流环的开环传递函数

2K2TTO(s) （2.3.6）

3R223R2R3s()s(2)s2TLTLTTL采用状态空间模型的能控标准型实现，传递函数（1）的状态空间表达式

AxBux yCxDu其中

（2.3.7）

3R23R2R2TL21TLTTL0C00A100B01002KD0 2T 5

2.3.2 系统的多胞型模型

考虑到在实际系统中的各种不确定因素，模型中的参数R, L在一定的范围内变化，引起系统矩阵A中第一行中的三个元素发生变化，因此系统（2）可以用一个多胞型模型表示如下

A(p)xBuxyC(p)x该系统的系统矩阵 S(p) （2.3.8）

A(p)B在以下给定的多胞型模型中取值 C(p)0

88S(p)Co1,S2S8akSk:ak0,ak1k1k1其中 S1A1C1BA2,S2C02BA8,,S8C08B是多胞型模型的8个顶点 0q1q2q3,q1q2r3,,r1r2r33R23R2Rq1r1q22 r2q32r3 （2.3.9）

TLTLTTL形成的8个系统矩阵。

2.3.3基于状态反馈的控制器设计

稳定三相电流型PWM整流器的直流侧电流是电流环的主要作用，所以我们考虑系统的两个性能，一是系统的瞬态响应应该有较小的震荡，二是系统的抗干扰能力。因此，这里考虑的指标一是闭环系统的极点，二是系统的H∞指标。

定理1[6]矩阵A的所有特征值均在半径为r，中心为(-q,0)的圆盘中的充分必要条件是存在对称正定矩阵X，使得线性矩阵不等式rXTqXXAqXAX0 rX1对于系统（1）的传递函数T(s)C(sIA)BD的H∞指标定义为

T(s)supmax(T(j))，既系统频率响应最大奇异值的峰值。

定理2[6] 对于系统（1），系统渐进稳定，且T(s)的充分必要条件是存在对称正定矩

6

ATPPAT阵P，使得线性矩阵不等式BBCPBIDCTDT0，结合以上分析，我们的目标是I对于系统（7）找到一个状态反馈控制律K，使相应闭环系统满足 1、 闭环系统的极点均在半径为r，中心为(-q,0)的圆盘中。 2、 闭环系统 T(s)尽可能小。

根据定理1,2，上述要求转化为带约束的凸优化问题

minrPstTqPPAP0对于系统多胞型模型，问题描述为

qPAP0rPATPPAPBCTTTBBID0CDI （2.3.10）

minqPAkP0rPATkPPAkTBkBkCkPBkIDkCTkDTk0I（2.3.11）

rPstTqPPAkP0k1,2,,82.3.4计算实例

为了验证以上方法的可行性，给出一个计算的例子。

在模型（7）中，设各参数的变化范围是 1To(s)在单位负反馈下的闭环传递函数

2 32s3.131s2.392s0.26112

s33.131s22.392s2.2611Tc(s)此闭环传递函数的单位阶跃响应见图2

7

1.4 1.2 System: sysC Peak amplitude: 1.18 Overshoot (%): 32.9 At time (sec): 3.92 System: sysC Settling Time (sec): 12.2 1 0.8 0.6 0.4 0.2 0 0

2 4 6 8 10 12 14 16 18 20

图2 校正前名义系统的闭环单位阶跃响应

可见，此时系统有较大的超调量，且调整时间很长。为了改善系统的性能，现在要求设计一个状态反馈律使得闭合极点在以（-3,0）为圆心，2为半径的圆盘内，且闭环系统的TC(s)尽可能小。根据参数R ，L的变化范围及（9）可以得到q1，,r1， q2，r2，, q3，,r3。再根据（10）得到一个有8个顶点的多胞型模型，对这个模型求解凸优化问题（11）得到状态反馈律K = [-5.4128 -3.6026 -0.5413]，保证在给定的参数变化范围内，系统的闭环极点均在在以（-3,0）为圆心，2为半径的圆盘内，且在参数的变化范围内，系统的TC(s)为了验证解的正确性，下面给出多胞型模型在8个顶点处的实际指标： 顶点1，闭环极点为{ -1.31 - 0.8439i ，-1.31 + 0.8439i -4.6377} , TC1(s)顶点2，闭环极点为{ -2.9407 - 1.2503i, -2.9407 +1.2503i, -4.4304} , TC1(s)顶点3，闭环极点为{ -1.8539 - 0.9523i, -1.8539 + 0.9523i, -4.9550} , TC1(s)顶点4，闭环极点为{ -2.3182 - 1.1206i, -2.3182 + 1.1206i -4.0265} , TC4(s)顶点5，闭环极点为{-2.2555 - 0.9155i, -2.2555 + 0.9155i, -3.9129} , TC1(s)顶点6，闭环极点为{-2.1384 - 1.0519i, -2.1384 + 1.0519i, -4.3861 } , TC1(s)顶点7，闭环极点为{ -2.9309 - 1.5968i, -2.9309 + 1.5968i, -2.5622} , TC1(s)顶点8，闭环极点为{ -2.7471 - 1.4365i -2.7471 + 1.4365i -3.1686} ，TC8(s)< 6.8415。

= 3.6945

= 3.6945 = 3.6945 = 3.6945 = 3.6945 = 3.6945 = 3.6945 = 3.6945

 8

根据多胞型模型的凸依赖性质，可知在所有的参数取值范围内，系统均满足要求。例如 顶点8处的归一化后的闭环传递函数TC8(s)应见图3

1.4 30.45，单位阶跃响32s8.663s27.02s30.451.2 System: sys Settling Time (sec): 2.05 1 0.8 0.6 0.4 0.2 0 0

0.5 1 1.5 2 2.5 3

校正后，系统无超调，且调整时间由校正前的12.2秒减少为2.05秒，可见系统性能得到明显的改善。

图3 校正后顶点8处的系统闭环单位阶跃响应

2.3.5 一些说明

稳定三相电流型PWM整流器的直流侧电流是电流环控制的的主要目的，所以我们考虑系统的两个性能，一是系统的瞬态响应应该有较小的震荡，二是系统的抗干扰能力。因此，这里考虑的指标一是闭环系统的极点，二是系统的H∞指标。对于这个系统来说，考虑到在实际系统中，参数L和R。由于对不同的顶点采用是相同的矩阵P，得到的多胞型系统的H∞指标的上限为6.8415，所以设计出的控制律有一定的保守性。正是这种保守性，提高了实际使用时的可靠性，因此对工程实际使用来说是有益的。实际计算得到，此多胞型系统在8个顶点处的TC(s)指标均为3.6945，实际上，对这个系统来说，参数L和R的变化不影响它

的H∞指标。进一步分析，可知系统H∞ 指标对应的频率是ω=0。所以，在实际工程设计中，只要考虑输入的直流分量即可。在本文中采用的设计方法的另外一个好处是，将一个系统的

9

指标设计转化为一个以线性矩阵不等式为约束条件的凸优化问题，所有可以很容易加入对系统其它指标的要求。

第二章 附录 线性矩阵不等式及求解方法

线性矩阵不等式是鲁棒控制理论中使用的主要数学工具之一，为了课程的完整性，这里专门作一介绍，本内容是作者根据有关参考资料整理而成。

A.1 线性矩阵不等式(LMI)的基本概念

线性矩阵不等式是指以下形式的矩阵不等式

F(x)F0xiFi0 (A.1.1)

i1m这里xRm是变量，x=(x1,,xm )，F0 =F0TR

nn

，Fi=FiTR

nn

(i=1,,m)是已知对称矩阵。这

(1)

(p)

时，称矩阵F(x)仿射(affinely)依赖于变量x。多个线性矩阵不等式F(x)>0,,F(x)>0等价于一个单一的线性矩阵不等式

diag(F(x)>0, , F(x))>0

(1)

(p)

(A.1.2)

因此从理论上来说，一个线性矩阵不等式与多个线性矩阵不等式无本质的区别。线性矩阵不等式的这一特性对本课题的研究有重要的意义，这也是我们选择它作为基本工具的主要原因之一。从下面的讨论中可以看到，只要将滤波器的单一指标的设计问题归结为一个线性矩阵不等式的求解问题，滤波器的多指标设计问题就可以归结为一个形如

diag(F(x)>0,, F(x))>0

(1)

(p)

的线性矩阵不等式的求解问题。

LMI (A.1.1)相当于关于x的m个多项式不等式。(A.1.1)是关于x的凸约束。也就是说，集合{ x | F(x)>0}是凸的。虽然从形式上看线性矩阵不等式(A.1.1)具有比较特殊的形式，实际上它可以表示一类很广泛的对变量x的凸约束。例如，凸线性不等式，二次不等式，矩阵模不等式以及在控制论中经常遇到的Lyapunov 和Riccati 不等式等等。 以下几类线性矩阵不等式在本文中经常用到，这里特别列出。由于本文不是专门研究线性矩阵不等式，这里仅不加证明地给出结论。 定理A.1.1（Schur补）线性矩阵不等式

T

T

Q(x)S(x)ST(x)R(x)0 (A.1.3)

其中Q(x)= Q(x) R(x)=R(x)，S(x)是形如(A.1.1)的矩阵，等价于非线性矩阵不等式

R(x)>0, Q(x)-S(x)R(x)S(x)>0 (A.1.4)

10

-1T

换句化说，非线性矩阵不等式(A.1.4)可以用线性矩阵不等式(A.1.3)表示。

根据以上定理A.1.1立即可以得到线性矩阵不等式

Z(x)I0 ZT(x)I等价于对矩阵模的约束‖Z(x)‖<1，其中Z(x)R

T

pq

。实际上，‖Z(x)‖<1等价于 I- Z(x)

Z(x)>0。这里，‖X‖=max(XTX)1/2,即矩阵的最大奇异值。

下面讨论以矩阵为变量的矩阵不等式，鲁棒控制理论中经常用到的两类矩阵不等式：

1、Lyapunov不等式

ATP+PA+Q<0 (A.1.5) 其中A,QR

nn

是已知常数矩阵，P=PT是变量。

3、 Riccati不等式

T

T

AP+PA+PBRBP+Q<0

T

T-1T

(A.1.6)

其中A,B,Q=Q,R=R>0是适维已知常数矩阵，P=P是变量。这是一个关于变量P的二次不等式，根据定理A.1.1它可以表示为线性矩阵不等式

ATPPAQPB0 (A.1.7) TBPR为了使上式符合一般的习惯和更有效地计算，不把矩阵不等式写作标准的(A.1.1)的形式，而

只是指出不等式中那一个是变量。例如，在不等式ATP+PA+Q<0中，如果变量是P，称为关于P的不等式。当然，矩阵不等式(A.1.5)可以归结为(A.1.1)所示的标准形式。下面看一个例子，考虑离散形式的Lyapunov不等式

ATXA-X+ Q<0，X=XTR2

这个矩阵变量的不等式可以写成(A.1.1) 的标准形式

ATXA-X+ Q=M0+x1M1+ x2M2+ x3M3

其中

M0=Q，Mi=ATEiA-Ei (i=1,2,3)

100100E1,E210,E301

00然而，(A.1.1)的标准形式从实用与计算的角度来看有不足之处。首先，按照(A.1.1)定义矩阵

序列M0,,Mn从计算的存储量来说是不经济的。当XRnn，每个Mi是nn矩阵，考虑到对称性，有n(n+1)/2个决定变量(decision variables)，因此需要n(n+1)/2个存储单元，总的存储单元与n4/4成正比，而存储相应的A，Q仅需要2n2个存储单元。其次，由于LMI具有不同的形式，变量X也有不同的结构，从标准形式 (A.1.1)转换为自然的形式 (A.1.3)的过程不能自动地完成，不利于计算机程序实现。因此，在实际计算中一般采用LMI的矩阵表达方式。

11

A.2 线性矩阵不等式标准问题

1、线性矩阵不等式问题（LMIP）

给出一个线性矩阵不等式F(x)>0，LMIP问题是找到一个可行解xf使得线性矩阵不等式F(xf)>0成立，或确定这样的xf不存在。 2、特征值问题（EVP）

极小化一个矩阵的最大特征值，满足线性矩阵不等式的约束，或确定在此约束下的解不存在，即

min  , s.t I-A(x)>0, B(x)>0 (A.2.1)

当A,B是变量x的对称矩阵时，约束是凸的，很明显，EVP问题等价于一个带有不等式约束的凸优化问题

min cx , s.t I-A(x)>0, B(x)>0

T

(A.2.2)

作为一个EVP问题的例子，考虑

ATPPAQPBmin  ,s.t P>0, 0 (A.2.3) TBPI其中ARnn, BRnP, CRmn, 是已知矩阵，是已知标量因子， ,P是待优化的变量。从上

面讨论可知，它等价于

min ，s.t ATP+PA+Q+-1PBBTP<0 (A.2.4)

3、广义特征值问题（GEVP）

广义特征值问题是极小化一对矩阵的广义最大特征值，满足线性矩阵不等式的约束，或确定在此约束下的解不存在，一般形式为

min  , s.t  B(x)-A(x)>0, B(x)>0，C(x)>0 (A.2.5)

这里A,B,C是仿射依赖于变量x的对称矩阵，上式也可以表示为

min max(A(x),B(x)) , s.t B(x)>0,C(x)>0 (A.2.6) 其中，max(X,Y)表示矩阵Y

-1/2

XY

-1/2

的最大特征值。GVEP问题是半凸(quasiconvex)优化问

题，因为约束条件是凸的，而目标函数max(A(x),B(x))是半凸的。也就是说对两个可行解x1，

x2以及标量0   1

max(A(x1+(1-) x2),B(x1+(1-) x2)) max{max(A(x1),B(x1)), max(A(x2),B(x2))}

(A.2.7)

4、 凸问题(CP)

min logdetA(x) ，s.t A(x)>0,B(x)>0 (A.2.8)

-1

-1

这里A,B是仿射依赖于变量x的对称矩阵，注意当A>0时，logdetA(x)是A的凸函数。作为一个CP问题的例子，考虑

min logdet P ，s.t P>0, viTPvi 1 (i=1, , L) (A.2.9)

-1

其中viRn为已知，P=PTRnn 是变量。下面我们给出这个问题的一个几何解释，用表示以原点为中心，由P确定的椭球 ={z|zPz  1}，约束条件是vi。由于椭球的体积正比于

T

12

det(P)

-1/2

,所以，极小化logdet P相当于极小化椭球的体积。求解(2.2.9)就可以得到以原点

-1

为中心，包含向量vi (i=1, , L)的最小椭球。即向量组vi (i=1, ,L) 的凸包(convex hull)。 无论理论上还是实际上，线性矩阵不等式的标准问题(LMIP, EVP, GEVP, CP)都可以有效地求解。也就是说，从理论上可以方便地证明这些问题解的存在性和求解算法的收敛性；在实际计算中有非常有效的算法进行求解。这里“求解”是指：确定这些问题的解是否存在，如果存在，找到一个可行解，它与全局最优解的误差小于预先给定的精度要求。

A.3 椭球算法(ellipsoid algorithm)

为了说明求解上述线性矩阵不等式的标准问题的思路，这里简单介绍一下椭球算法。这个算法给出了求解上述问题的一个简单而明确的思路。当然，从实际计算的角度来看，近几年发展起来的内点算法(interior-point algorithm)更加有效，该算法我们将在后面介绍。

椭球算法基本思想如下，假设要求解的问题至少有一个最优点，在求可行解时，认为可行点是一个最优点。从一个椭球0开始，假设0中至少有一个最优点。计算一个切平面(cutting plane)，这个平面经过椭球0的中心x0，也就是说，找一个非零向量g0使最优点落在半空间{z| g0(z- x0)  0}中（后面我们解释对每一个标准问题，如何找向量g0）。从而，被分割的半椭球 0 {z| g0(z- x0)  0}至少包含一个最优点。接下来，计算包含这个半椭球的最小的椭球1，显然，1中至少包含一个最优点。继续这一过程，可得到一系列的体积越来越小的椭球，当椭球的体积小于预先给定的精度要求时，从最后得到的椭球中任取一点就是要求的最优解。

下面给出具体的计算公式。一个椭球可以表示为

其中A=A>0。包含半椭球

的体积最小的椭球是

这里

1={z| (z-a1)A1 (z-a1)1}

T

-1

T

T

T

={z| (z-a)A (z-a)1}

T-1

(A.3.1)

={z| (z-a)A (z-a)1, g0(z- x0)0}

T-1T

(A.3.2)

(A.3.3)

Ag1m22Ta1a,A12(AAg1g1A),g1m1m1m1ggAgT

此公式当m2时成立。当只有一个变量时，包含半区间的最小区间就是它本身。这时椭球

算法退化为我们熟悉的对分法。

椭球算法以初值x,A开始，按下述过程进行 计算过x(k)的切平面g(k)；

(A.3.4)

这个过程产生一系列的椭球，并保证每个椭球中都包含一个最优点。它们的体积按几何级数

13

(0)

(0)

递减 vol(

(k)

) e

-k/2m

vol(

(0)

)，这里vol( )表示椭球的体积。

下面讨论标准问题中切平面的计算方法。 1、LMIP问题 考虑LMI

F(x)F0xiFi0

i1m如果x不满足这个不等式，存在一个非零的u使

T

uF(x)uu(F0xiFi)u0

TTi1m (A.3.5)

令gi=-uFiu，对任意z满足g(z-x)0有

T

uF(z)uu(F0xiFi)uuTF(z)ugT(zx)0

TTi1m因此，每一个可行解均在半空间{ z | g(z-x)< 0}内，所以，g确定了一个LMIP在点x的切平面。

2、EVP问题 考虑EVP

min cx s.t F(x)>0

T

T

(2.3.6)

假设点x使 F(x)<0,我们可以用以上在LMIP中的方法构造一个切平面。这时，我们去掉半空间{z|g(z-x)0}，因为这里面的点都不满足不等式 F(x)>0。如果给出的点x满足 F(x)>0，这时，g=c确定了一个LMIP在点x的切平面。去掉半空间{z|g(z-x)>0}，因为这里面的点(无论是不是可行解)的目标函数的值均大于x，因此不可能是最优解。 3、GEVP问题 考虑以下公式

min max(A(x),B(x)) ,s.t B(x)>0, C(x)>0 设所给出的x是可行解，取任意向量u0满足

(max(A(x),B(x)) B(x)- A(x))u=0 由下式定义的g

g=u(max(A(x),B(x)) B(x)- A(x))u

即为一个GEVP在点x的切平面。因为，对任意z

u(max(A(x),B(x)) B(x)- A(x))u=g(z-x)

因此，如果g(z-x)0，则有max(A(z),B(z)) max(A(x),B(x))。

T

T

T

T

T

T

(A.3.7)

任给一个点x，如果不是可行解，可以用以上在LMIP中的方法构造一个切平面。现在，假

14

4、CP问题

最后考虑CP问题，任给一个点x，如果不是可行解，可以用以上在LMIP中的方法构造一个切平面。现在，假设所给出的x是可行解，切平面由目标函数的梯度给出

-1

logdetP(x)logdet(P0xiPi)

i1m (A.3.8)

即 gi=Tr(PiP(x))。由于目标函数是凸的，对所有的z

logdetP(z)1logdetP(x)1gT(zx)

特别地，g(z-x)>0 意味着

所以z不可能是最优解。

T

logdetP(z)1logdetP(x)1

(A.3.9)

A.4 解正定规划的内点方法(interior-point method)

从1988年以来，内点方法被用来求解线性矩阵不等式的标准问题。目前，该方法已经非常成熟，经验证是一种高效率的求解线性矩阵不等式的标准问题的方法。本节中我们简要介绍基本概念。

A.4.1 线性矩阵不等式的的解析中心(analytic center)

线性矩阵不等式解析中心的概念在内点方法中起着重要的作用，这里简述一下它的概念，考虑线性矩阵不等式

F(x)F0xiFi0

i1m这里xR是变量, F0=F0T, Fi=FiT R的可行解集合

mnn

(i=1,,m) 是已知的对称矩阵。以X表示该不等式

X={ xRm | F(x)>0 }

显然，函数

(A.4.1)

仅当xX时有确定的值。当x靠近X的边界时，该函数趋于无穷大，也就是说，它是X的

障碍函数(barrier function)。

假设X是非空的有界集合，这意味着矩阵 Fi=FiR

T

nn

, (i=1, ,m)线性无关(否则，X将包

(A.4.2)

含一条直线)。可以证明，(x)是X的严格凸函数，所以它有一个唯一的最小点

x*=argmin(x)

称x*为线性矩阵不等式F(x)>0的解析中心。等价地有

x*=argmax detF(x)

15

也就是说，在所有的形如F(x)的正定矩阵中，F(x*)具有最大的行列式。

下面讨论线性矩阵不等式的解析中心计算方法，实际上，这是一个特殊的CP问题。只要选择合适的步长，牛顿迭代法可以有效地计算线性矩阵不等式的解析中心x*。在X中任取一点，采用迭代公式

(k)

x

(k+1)

=x-aH(x)g(x)

(k)

(k)

(k)(k)(k)-1(k)

(k)

(A.4.3)

这里，a为第k次迭代的阻尼因子，H(x)和g(x)分别代表(x)在x处的海塞(Hessian)阵和梯度。阻尼因子的取法如下

(k)

111(x(k))(x(k))1/4(x(k))1/4

(A.4.4)

(x(k))g(x(k))TH(x(k))1g(x(k))可以证明，这样选择的阻尼因子总会使x

(k+1)

X，并且保证x收敛于x*。

(k)

A.4.2 正定规划(PDP)问题与LMI

考虑

min cx , s.t F(x) 0

T

(A.4.5)

这里F(x)是x的仿射函数，(2.4.5)称为一个正定规划(PDP)问题，当约束是凸的时，称(A.4.5)是一个凸规划。对xR如果F(x) 0，称x为可行解，如果F(x) >0，称x为严格可行解。 记PDP(A.4.5)的最优值为r

optoptm

r=inf{ cx | s.t F(x) 0}

optT

opt

(A.4.6)

当PDP(A.4.5)无下界时，r = -。习惯上，当PDP(A.4.5)无解时，我们认为r =+。求解PDP(A.4.5)的含意是

(1) 如果该问题无解(r =+)，给出无解的证明。

(2) 如果r 为有限值，给出一个可行解x满足r cx-，称x为-次优解，这里>0。 (3) 如果目标函数无下界，即r =-，给出一个可行解x和一个可行方向v，沿着方向v目标函数趋于-。

下面讨论LMI与PDP的关系，首先多个LMI可由块对角矩阵合成一个大的LMI。G1(x)  0,,GL(x)  0等价于F(x)=diag(G1(x)  0, ,GL(x) )  0,所以，PDP问题

opt

opt

opt

T

opt

min cx , s.t Gi(x)  0 ( i=1, ,L)

T

(A.4.7)

16

等价于PDP(2.4.1)。考虑凸优化问题

min f(x) , xC (A.4.8)

其中f：RkR是凸函数，C是凸集。 如果可以找到一个仿射函数Fo：Rk+1Rpp使

f(x)  t  Fo(x,t)  0 (A.4.9)

kqq

和仿射函数Fc：RR使

x C  Fc(x)  0 (A.4.10)

则凸优化问题 (2.4.8)可以表示为一个PDP

min t , s.t Fo(x,t)  0, Fc(x)  0

在这个PDP中，m=k+1,xR ,tR，向量c1=, ,=ck=0,ck+1=1, F(x) R

k

nn

，F(x)  0，n=p+q，

F(x)=diag(Fo(x,t), Fc(x))。

LMI(A.4.9),(A.4.10)称为函数f及约束C的正定表示(PDR)。为了将一个凸优化问题转化为一个PDP必须找到它的目标函数与约束的PDR。下面看几个例子。 (1) 线性约束 A  xb，AR

nn

，bR，(向量不等式)。相应的PDR是

T

n

F(x)  0，F(x)=diag(Ax-b)

因此，它可以表示为PDP

min cx , s.t Ax  b

T

(2) 凸二次函数

任意凸二次函数可以表示为f(x)=(Ax-b)(Ax-b)-cx-d,它的PDR是

T

T

IF(x,t) (Axb)TAxb0

tcTxd(3) 欧氏 (Euclidean) 模

函数‖Ax-b‖的PDR是

tIF(x,t)T(Axb) pq

Axb0 t(4) 矩阵的模 矩阵A(x)  R

，A(x)是x的仿射函数。‖A(x)‖F的PDR是

tIF(x,t)TA(x)(5) 矩阵的最大特征值

T

A(x)0 tI矩阵A(x) =A(x)是x的仿射函数。最大特征值 max{A(x)}的PDR是

F(x,t)=tI-A(x) 0

(6) 最大值与交集

设F1(x,t)  0, F2(x,t)  0是凸函数，分别是f1 , f2的 PDR，则f=max(f1,f2)的PDR是

17

diag(F1(x,t), F2(x,t))  0

类似地，设凸集C1,C2的PDR是F1(x)  0, F2(x)  0，交集C1C2的PDR是

diag(F1(x), F2(x))  0

(7) 函数和集合的和(sum)

设F1(x,t)  0, F2(x,t)  0分别是f1 , f2的PDR，则f= f1+f2的PDR是

F1(x,v)  0, F2(x,t-v))  0

其中，v是一个新变量。类似地，设凸集C1,C2的PDR是F1(x)  0, F2(x)  0，和集C1+C2的PDR可以表示为

F1(w)  0, F2(x-w))  0

其中，w是一个新变量。

A.4.3 对偶问题

定义A.4.1 PDP min cx , s.t F(x)  0的对偶问题为

max -TrF0Z s.t TrFiZ=ci, Z  0 (A.4.11)

对偶问题也是一个PDP问题，可以表示为原始问题的方式。事实上，如果F1, ,Fm线性

{Z | Z=Z R

可以表示为

{G(y)=G0+ y1G1+, ,+ ypGp |y  R }

p=n(n+1)/2，G i为适维矩阵，设d R，d= TrF0Gi 则有

cy= TrF0(G(y)- G0)

对偶问题变为

min dy , s.t G(y)  0

T

Tp

p

T

nn

T

，TrFiZ=ci}

(A.4.12)

它具有与原始问题相同的形式。

对偶PDP的一个关键特性是它给出了原始PDP问题的最优解的一个下界，对所有的可行解x和Z有- TrF0Z  cx，这是我们下面用到的对偶内点方法的关键。由于x和Z是任意的，所以对原始问题的最优解r

opt

T

也有-TrF0Z  r。在解决实际问题时，大多数算法是产生一个

opt

opt

最优解的近似值xS，保证c xS- r  ，为预先给定的精度要求。 定义A.4.2如果x是一个原始问题的可行解，Z是对偶问题的可行解，称

= cx+ TrF0Z

T

T

(A.4.13)

为关于可行解x和Z的对偶区间，显然有

-TrF0Z  r  cx

18

opt

T

是包含最优解r定理 A.4.1 设

opt

的区间。下面我们以定理的形式给出对偶问题的解的一些性质。

p*=infimum cx s.t x  R ,F(x)  0 d*=supremum -TrF0Z s.t TrF0Z ,Z  0

则p*=d* ，只要下列四个条件之一成立 (1) 原始问题有严格可行解 (2) 对偶问题有严格可行解

(3) 原始问题的可行域X={x|F(x)  0, cx =p*}非空且有界

(4) 对偶问题的可行域Z={Z | TrFiZ=ci, (i=1,,m), Z  0, -TrF0Z =p*}非空且有界。 定理 A.4.2 设向量xR ,F(x)  0是原始问题的最优解，如果存在一个矩阵Z  0，且TrFiZ=ci (i=1, ,m),cx=-TrF0Z，则有以下结论：如果cx+TrF0Z  ，则xR ,F(x)  0是原始问题的次优解。

反之，设矩阵Z 0，且TrFiZ=ci , (i=1,,m)是对偶问题的最优解，如果存在一个向量xR 满足F(x)  0，cx=-TrF0Z，有以下结论：如果cx+TrF0Z  ，则Z是对偶问题的次优解。 定理 A.4.3 集合{ xR ,F(x)>0}是空集，当且仅当存在

Z  0，Z  0，TrFiZ=0, (i=1, ,m)，TrF0Z  0

m

m

T

T

m

T

T

m

m

T

Tm

集合{ xR ,F(x)  0}是空集，当且仅当存在

Z  0，Z  0，TrFiZ=0, (i=1, ,m)，TrF0Z < 0。



定理A.4.4 存在向量v满足

vFii1mi0,cTv0当且仅当集合

{ Z | Z>0, TrFiZ=ci (i=1, ,m) }

是空集。类似地，存在向量v满足

vFii1mi0,cTv0当且仅当集合

{ Z | Z0, TrFiZ=ci , (i=1, ,m) }是空集。

A.4.4 原始—对偶方法

这个方法的基本思路是求解 PDP

min cx+TrF0Z

s.t F(x)  0 ,Z  0, TrFiZ=ci , (i=1, ,m)

(A.4.14)

显然，这个问题的最优值为零。下面我们给出求解这个PDP的公式。定义位能函数(potential function)

19

T

1/2

(x,Z)=qlog(cx+TrF0Z)-logdetF(x)-logZ-nlogn

T

(A.4.15)

其中，q=n+vn，参数v是式中第一项的权重，将 (x,Z)分解为两项

(x,Z)= (vn)log(cx+TrF0Z)+ (x,Z)

1/2

T

(A.4.16)

(A.4.16)中的第一项是目标函数值减少的度量，第二项表示变量的当前值与中心的偏差。也就是说，(A.4.16)中的第一项只依赖于对偶区间。当对偶区间0时，

(vn)log(cx+TrF0Z)  -

因此，它的减少相当于对偶区间的减少。第二项(x,Z)表示变量的当前值与中心的偏差，当(x,Z)趋于正定锥(positive-definite cone)的边界时，(x,Z)  +。因此，它的减少相当于(x,Z)向解析中心靠近。解 (A.4.14)相当于求(2.4.7)的最小值。求解的方法可分为两类，一类是利用改进的牛顿迭代法，另一类方法是将问题转化为最小二乘问题。原始—对偶方法的基本思想是产生原始—对偶变量的迭代满足

(x(k+1), Z(k+1)) (x(k), Z(k))- , >0

注意，函数 (x, Z)中的第一项是凹的，为了能够利用凸优化技术，将函数 (x, Z)中的第一项在x(k), Z(k)处线性化，得到一个凸函数c(x, Z)，这个函数具有以下特性

c(x, Z)   (x, Z), c(x(k), Z(k))   ( x(k), Z(k))

因此，如果选择x(k+1), Z(k+1)使c(x, Z)减少某个数量，(x, Z)减少的更多。

下面的方法实际上是一种多变量的拟牛顿迭代方法，在牛顿迭代公式中使用(x,Z)的梯度，但是在计算(x, Z)的海塞(Hessian)阵时是近似的，即不考虑(x, Z)中的凹项。下面给出计算公式。将(x, Z)的第一项线性化得到 c(x, Z)= ( x(k), Z(k))+p(cT(x- x(k))+TrFo(Z -Z(k))-logdetF(x)-logdetZ-nlogn (A.4.17)

T(k)(k)

其中 p=q/( cx+TrFoZ)，注意在c(x, Z)中，x, Z是分离的。将c(x, Z)按x, Z写成两项

其中

c p(x) = pcTx(k)-logdetF(x) c d(Z)= pTrFoZ-logdet Z

因此，在计算牛顿方向的增量xp ,Zd时，原始问题与对偶问题可以分开考虑。注意到

mm1m111(xv)(x)TrF(x)FiviTrF(x)FiviF(x)FjvjFi12i1j111/2T

c(x, Z)= cp(x, Z) +cd(x, Z)+const (A.4.18)

(x)viTrF(x)1Fii1m1vivjTrF(x)1FiF(x)1Fj2i1j1

mm

(A.4.19)

(x)的梯阵度g(x)和海塞阵H(x)分别为

gi(x) =-TrF(x)-1Fi = -TrF(x)-1/2Fi F(x)-1/2 (i=1,,m)

Hi j(x) = TrF(x)-1Fi F(x)-1 Fj = Tr( F(x)-1/2Fi F(x)-1/2)( F(x)-1/2Fj F(x)-1/2 )

(i , j = 1,,m) 因此有

20

m1m1mxpargmincvTrFFiviFiviFFjvjF1i12i1j1vRmT1argmincTvviTrF1FiTrF1vRmi1m1vivjFiF1FjF12i1j1mm

 (A.4.20)

令Gi=F-1/2Fi F-1/2 (i=1,,m)，CRnn是满足 TrFiC=ci ，(i=1,,m)的任意矩阵（例如可取C=Z），则xp可以简化为

(A.4.21)

这是具有m个未知数n(n+1)/2个方程的最小二乘问题。类似地，在所有的可行方向中，Zd 是使二次近似cd(Z)最小的方向。

1dZdargminpTrF0VTrZ1VTrZ1VZ1

2V:TrFiV0i1,,m (A.4.22)

这是一个最小二乘问题，引入Lagrange乘子[10]，vi (i=1,,m) (2.4.22)可以表示为

F0-Z-1+Z-1Zd Z-1+( v1F1+ +vmFm)=0 Tr FiZd=0 (i=1,,m)

或等价地，记Gi=Z1/2Fi Z1/2 (i=1,,m)

Z1/2F0 Z1/2-I+Z-1/2Zd Z-1/2+( v1G1++ vmGm)=0 (A.4.23) Tr Gi Z-1/2Zd Z-1/2=0 (i=1,,m)

所以，vi (i=1,,m)是以下最小二乘问题的解

minZF0ZmvR1/21/2IGivi

i1Fm (A.4.24)

这里‖‖F表示Frobenius模，XFX1/2表示矩阵X=XT>0的正定平方根因子。(TrX2)1/2，

(A.4.23)中有n(n+1)/2个等式，m个变量。解出vi (i=1,,m)代入(A.4.23)可得到Zd。综上所

述，得到基本的原始—对偶算法(BPDA) 给出严格的可行解x, Z（初始点），重复以下步骤 (1) 解最小二乘问题(A.4.21)求xp

(2) 解最小二乘问题(A.4.24),由(A.4.23)求Zd (3) 求p,qR，min ( x+xp, Z+Zd) (4) 更新 xx+xp, ZZ+Zd 直到满足 cTx+TrFZ<。

下面讨论这个算法的收敛速度。记

pFi1m1/2FiF1/2xpi,dFZi1m1/2ZdZ1/2

F (A.4.25)

则有以下结果。

定理A.4.5 设x, Z是严格的可行解，xp，Zd是相应的牛顿方向，记 p=1/(1+p), q=1/(1+d) 则 ( x+pxp, Z) ( x, Z)- p+log(1+p)

21

( x, Z+qZd)  ( x, Z)- d+log(1+d)  定理A.4.6 max{p ,d }0.35。

从以上两个定理中不难看出，每一次迭代至少可使位能函数减少

log0.35-log(1-0.35)=0.78 (A.4.26)

所以，该算法的收敛速度是O(n1/2)。

可以对上述BPDA算法进行改进，使每次迭代只须求解一次最小二乘问题。基本思想是考虑(A.4.21)式，xp使(A.4.21)取最小值的条件是

pciTrFFixpjTr(FiF1FjF1)0

1j1m也就是说，矩阵

(1/p)TrFi(FxpjTr(FiF1FjF1)ci

1j1m满足对偶可行解的等式约束，所以

Zp(1/p)(FdxpjTr(F1FjF1)Z

1j1m(A.4.27)

是一个对偶可行方向。改进后的算法以这个Zp为对偶搜索方向，在以xp,Zp定义的平面上搜索最优点，具体算法如下(NNPDA) 给出严格的可行解x, Z（初始点），重复以下步骤 (1) 解最小二乘问题(2.4.21)求xp (2) 由(2.4.27)求Zp

(3) 求p,qR，min ( x+xp, Z+Zd) (4) 更新 xx+xp, ZZ+Zd 直到满足 cTx+TrFZ<。

综上所述，对偶--原始方法的基本思想是：首先求得位能函数的可行方向xp, Zd（牛顿迭代法或最小二乘法），然后在可行方向定义的平面上搜索得到目标函数的最优点，重复这一过程直到对偶区间小于给定的精度要求。

22