优化设计-2

第1章 绪论

1.1 最优化问题的提出

“优化”既是一个专业术语，也是一个通俗的词语，这一方面说明优化问题的广泛性；另一方面也说明解决优化问题具有一定的难度，需要有专门的理论和技巧。优化问题来源于求某一设计（广义的设计）的最优结果，用数学观点来说就是求用某一指标或某几个指标描述的设计的最大值或最小值。设计的决策包含优化的过程，其中有通过以往经验判断得出的决策，有通过枚举或多方案比较得出的决策，而经济的做法则是通过对设计建立数学模型，通过解析或数值计算寻找到决策的依据，用以指导设计的实施。例如，某设计的模型可用一元函数f（x）来表示，对其进行最优化设计就是求该一元函数的最大值和最小值。如果一元函数是单调函数，则函数的最大值或最小值会在变量x的边界上取得；如果一元函数是二次多项式，则函数的极值在函数曲线的顶点上取得；如果一元函数是高次多项式，函数曲线有多个极值点，则求函数的最大值或最小值问题就变得复杂起来。对多元函数的极值问题更是如此，需要用到局部或全局优化算法来求解。

线性规划问题是目标函数和限制条件都是线性函数的问题，在生产和生活中很多问题都可抽象为线性规划问题。下面以线性规划举例说明优化设计问题的提出、建模及求解的全过程。 例1

有一名学生，期末有5门功课要考试，可用的复习时间有18h。假定这5门功课分别为数学、英语、计算机基础、画法几何和专业概论。如果不复习直接参加考试，这5门功课预期的考试成绩分别为65分、60分、70分、60分和65分。复习以1h为一单位，每增加1h复习时间，各门功课考试成绩就有可能提高，每复习1h各门功课考试成绩提高的分数分别为3分、4分、5分、4分和6分。问：如何安排各门功课的复习时间可使平均成绩不低于80分，并且数学和英语成绩分别不低于70分和75分？ 解：

设分配在数学、英语、计算机基础、画法几何和专业概论这5门功课的复习时间分别为x1，x2，x3，x4，x5，则可列出如下的目标函数和限制条件： min f（x）=x1+x2+x3+x4+x5 s.t. x1+x2+x3+x4+x5≤18

(3x1+4x2+5x3+4x4+6x5+320)/5≥80 3x1+65≥70 4x2+60≥75 3x1≤35 4x2≤40 5x3≤30 4x4≤40 6x5≤35

x1,x2,x3,x4.x5≥0

根据所给列出的主要方程。但是根据实际情况，各门功课的成绩不能大于100分，各门功课的复习时间不能是负数，因此还需补充这几个限制条件。由此看出，这是一个在满足限制条件（约束条件）的情况下，求最少复习时间的问题。下面用MATLAB优化工具箱求线性规划的函数linprog（）来求解此问题。 线性规划问题数学模型

min CTX s.t. AX≤B AeqX=Beq Lb≤X≤Ub

其中：C，X，B，Beq，Lb，Ub为向量；A，Aeq为矩阵。 例2

资源分配问题是线性规划中的一类问题。这里所说的资料其含义广泛，可以使一般的物质资料，也可以是人力资源。资源分配问题可描述为生产若干种产品（广义的产品）需要几种不同的资源，如原料消耗量、设备使用量、人力需要量等。各种资源供应量有一定限制，所生产的产品有不同的利润或花费不同的费用。所求问题是在所消耗资源和资源供给量限制的条件下，求生产不同的产品的数量，使收益最大或费用最低。

某工厂要生产两种规格的电冰箱，分别用Ⅰ和Ⅱ表示。生产电冰箱需要两种原材料A和B，另需设备C，生产两种电冰箱所需原材料、设备台数、资源供给量及两种产品可获得的利润如表所示。问：工厂应分别生产Ⅰ，Ⅱ型电冰箱多少台，才能使工厂获得最多？

资源 设备C 原材料A 原材料B 单位产品获利 电冰箱Ⅰ用原材料限制 Ⅰ 1 2 0 220 Ⅱ 1 1 1 250 ≤800kg 资源限制 1200台时 1800kg 1000kg 求最大收益 解：

设生产Ⅰ，Ⅱ型电冰箱的数量分别是x1，x2，则可获得的最大收益为 max f(X)=220x1+250x2, X∈R2 s.t. x1+x2≤1200

2x1+x2≤1800 x2≤1000 x1≤800 x1,x2≥0

通过以上两个例子，我们初步了解了优化设计求解的过程，以及优化设计的“威力”。 例1-2中使用了标准化得优化数学模型，而优化问题数学模型的一般描述为 min(max) f(X)=f(x1,x2,x3,…,xn), X∈Rn s.t. gu(x) ≤(≥)0, u=1,2,…,l hv(x)=0,v=1,2,…,m

目标函数和约束条件可以是线性的，也可以、是非线性的。

1.2 最优化问题分类

工程实际问题多种多样，依据不同的分类方法，它们属于不同的类型，相应地有不同的解法。例如根据问题的性质，可将问题分为静态问题、动态问题、确定性问题、随机性问题、模糊问题、连续问题、离散问题、逻辑问题等。求解问题时首先要根据问题遵循的基本定律建立相应的数学模型，模型的类型与问题的类型往往是一致的。不论是那一种类型的问题或模型，其求解可通过实验法、解析法或数值方法来实现。随着计算机技术及各种数值方法的发展，利用计算机求解实际问题越来越普遍和快捷。常遇到的工程问题其数学模型一种是代数模型；另一种是用常微分或偏微分方程表示的模型，如弹性体静力平衡微分方程。有限元法是求解偏微分方程非常有效的方法，其理论基础就是求能量函数或范函的极值。

对于优化设计问题来说，若目标函数和约束函数都是线性函数，则这样的问题就是线性规划问题，，如例1-1和1-2中数学模型；若目标函数或约束函数中含有非线性函数，则这样的问题是非线性优化问题。

有些优化问题有约束条件，而有些优化问题没有约束条件，据此可以将最优化问题分为无约束问题和有约束问题。无约束问题在经典优化设计中占有重要地位，其求解方法是某些约束优化问题求解的基础。

此外，目标函数中变量（称为设计变量）可能只含有一个，也可能含有多个，相应的优点问题就分别称为单变量问题和多变量问题。单变量优化问题的解决称为一维搜索方法。经典多变量优化算法中确定搜索方向最优步长的问题就是一维搜索问题。因此，一维搜索算法是多变量优化算法的基础。

根据目标函数的多少，最优化问题又可以分为单目标函数问题和多目标函数问题；根据设计变量取值的性质，最优化问题又可分为单目标函数问题和多目标函数问题；根据设计变量取值是否连续，最优化问题又可以分为连续优化问题和离散优化问题。

根据求解算法是否含有导数运算，可以将优化算法分为含导数的优化算法和不含导数的优化算法，不含导数的优化算法又称为直接算法。

离散优化问题通常称为规划问题，如资源配置、生产管理、最短邮路等问题。一些启发式优化算法如遗传算法、粒子群算法不仅适合于求解连续优化问题，也适于求解离散优化问题。下式是最优化问题数学模型：

这也是MATLAB优化工具箱函数中优化数学模型采用的形式。

经典优化算法大多属于线搜索方法，即从某一初始点出发x(0)，沿搜索方向d(0)，按一定的步长α(k)在约束条件限定的范围内进行搜索，一般的迭代格式为

x(k+1)=x(k)+α(k)d(k)

根据搜索方向d构造方法的不同，就形成了不同的优化算法。

1.3 优化模型的图形表示

由图1.1可知，设计变量的取值被限定在约束函数规定的区域内，满足约束条件的最优解位于约束函数曲线的交点上。将优化模型用图形表示就是优化问题的图解法。MATLAB有方便、灵活的绘图函数，对一些二维或三维优化问题应用绘图函数可以帮助了解目标函数性状和约束空间。对一些实际问题，通过图形表示可以了解设计变量的取值范围，也可以通过图解直接得出优化问题的解。先介绍MATLAB常用的绘图函数，然后通过典型优化问题的分析了解基于MATLAB的图解法。

MATLAB绘图函数包括平面曲线函数plot()、三维曲线绘图函数plot3()、三维曲线绘图函数mesh()和surf()，以及将三维曲面投影到平面上的取等值线函数contour()。要绘制较完美的图形，还要对曲线线型、线宽、线的颜色、绘图点标记的形状、标记边框颜色、标记填充颜色等进行定义，此外还要对坐标轴刻线、坐标轴名称、坐标轴取值范围、曲线图例等进行说明，对所绘图形修饰性的说明或定义也可在图形绘制完成后在图形显示窗口通过编辑命令来完成。

1. plot()函数

是在直角坐标系中绘制平面曲线的基本函数，要求输入的参数是横坐标x和纵坐标y的值。x和y用行向量或列向量来表示。 例1-3 绘制下面函数曲线： y(x)=2sinx+lnx

用到的有关函数有内置函数定义函数inline() 、坐标点生成函数linspace()、函数取值函数feval()。坐标点生成函数linspace()用于生成等距点，该函数有3个输入参数，前两个参数说

明坐标的起止点，第3个参数表示生成的离散点数，总的点数包括两个端点。如linspace(1,2,3)，结果为[1.0 1.5 2.0]，即区间[1 2]上生成3个点，在给定区间内插入一个点1.5.若不指定第3个参数，则自动认为生成100个点。生成等距坐标点的另一种方法是用分割符来实现，如x=1：0.1：10，结果是从1开始，按增加0.1递增，直到小于或等于终点。除可生成等距点外，还可以生成对数分隔点，所用函数为logspce()，该函数也有3个输入参数，其中前两个参数是以10为底的指数，第3个参数与linspace()函数相同。绘图时，函数内向量元素的个数要相等，计算向量长度、矩阵维数的函数分别是length()和size()。

二维绘图函数还有简洁绘图函数ezplot、绘直方图函数bar、在图形加上误差带数errorbar、用给定精度绘图的函数fplot、绘极坐标图函数polar、绘统计分布图hist、绘台阶图函数stairs、绘极坐标统计分布图函数rose、绘火柴杆图函数stem、绘图区填充函数fill、绘矢量场图函数quiver。 2. plot3()

用来绘制三维曲线，它有3个输入参数，分别对应直角坐标系中的x，y，z方向的坐标点，各坐标点以向量表示。 3. mesh()

用来绘制网格状三维曲面，它有3个输入参数，每个参数都是二维矩阵，它们分别是与平面矩形区域网格点对应的x，y，z坐标值。因此，要绘制三维曲面，首先要在xy平面上根据x，y的取值范围划分网格，然后再计算网格点上与x，y对应的z坐标值。网格坐标生产函数为meshgrid()，该函数根据x，y方向的两个一维向量生成平面网格。 4. contour()

三维曲面等值线函数，该函数常用格式：

contour(z) 根据矩阵z绘制三维曲面x，y坐标用矩阵行、列的序号表示。 contour(z,n) 绘制n条等值线。

contour(z,v) 按增量v绘制等值线，v为向量。

contour(x,y,z) 用给定的x，y坐标在xy平面上绘制z等值线。 contour(x,y,z,n),contour(x,y,z,v) 其中参数n，v含义与上同。

contour()函数既可以在三维图形中，也可在二维图形中。应用该函数结合优化问题的约束条件可进行优化问题的图解分析。

例1-4 用图形表示如下优化模型，并求解： min f(x)=4x1^2-x2^2-12 s.t. x1^2+x2^2=25 (x1-5)^2+(x2-5)^2-16 解：

绘制目标函数曲面、约束函数曲线及求解程序如下。

从图1.3可以看出，约束函数曲线在点（1，5）处相交，该点就是目标函数满足约束条件的解。应用MATLAB优化工具箱函数fmincon()得出的结果为X*=[1.0013 4.8987]T。通过图解分析可对优化问题的几何性质有更清晰地认识，从而可得出更为精确的解。

三维绘图函数中surf函数和riobbon函数也很有用，它们分别用来绘制带阴影的三维曲面以及用条带表示的三维曲面或曲线。

除可利用MATLAB进行函数曲线绘制和分析外，其他的一些数学软件也能轻松完成类似的工作，其中应用较广的是Maple，MathCAD和Mathematics。

1.4 有限元法引例

有限元法是求解偏微分方程的一种数值方法，在解决弹性力学、流体力学等学科的问题中得

到广泛应用，很多结构或流体分析商业软件都采用有限元法求解。下面通过阶梯轴有限元列式的建立来说明有限元法的基本思想和基本格式的推导过程。 一受轴向力作用的阶梯轴，采用虚位移原理建立该阶梯轴的有限元列式。虚位移原理可叙述为：弹性体在给定的约束和一组外力作用下处于平衡时，外力在约束和变形协调条件允许的任意一组虚位移上所做的虚功等于弹性体内应力在相应虚应变上所做的功，即

?*TF???*T?dv

v对阶梯轴的典型单元，虚位移原理可表示为Fi?ui?Fj?uj??le0?*T?Aedx

du dx应力与应变的关系为 （本构关系）??E?

位移与应变的关系为 （几何关系）??设位移模型为u??1??2x 对第①单元，边界条件为x?0将边界条件代入位移模式，得

u?uix?leu?uj

?1?uie?2?uj?uilexxu?ui(1?)?uj?(Nilele?ui?Nj)??

?uj?记????ui???ui?*??，则虚位移???。

u?j???uj?dNj??ui?1??(???dx?le??uj?1?ui?)???B?e le?uj?du?dNi??根据几何关系，有??dx??dx虚应变为??B?

**?*T??*TBT?(?ui?1????le???uj)??

1????le??1?ui?)?? le?uj?*T根据本构关系，有

??E??EB?e?E(?1le将以上各式代入虚功方程，得?F??le0?*T?Aedx

1?ui?)??Aedx le?uj?即(?ui?uj)????(?ui0F?j??Fi?ele?1????1?l??uj)?e?E(?1le????le??

因为δ*是任意一组可能的虚位移，要使上式相等，必有

le?Fi?????0?Fj?e?1????1?le?E(??1?le????le??1?ui?)??Aedx le?uj?即

?1e?2?Fi?le???E??1?Fj???l2?e1??1?2??2le?le?ui?lAe???dx?EAe?e1?0?uj??1?22??le??le1??2?le?1?le2??e?ui????uj?eFe?ke?e,?EAe?2lek??e?EAe??l2e?(1)k12(1)k22?EAe??le2? EAe?le2??下面进行总体刚度矩阵的装配。扩展单元刚度矩阵为

K~~(1)(1)?k11?(1)??k21?0?00??0?,0??~~(1)K?K~~(2)?00(2)??0k11?(2)?0k21?0?(2)? k12?(2)?k22?总体刚度阵为

K?K~~(2)

则整体有限元列式为KU?F,U?[u1u2u3]T,F?[F1F2F3]T

1.5 多学科设计优化集成软件iSIGHT简介

是一个集成优化平台，它一方面将各学科的专业知识和分析软件集成，另一方面提供丰富的优化探索策略，包括各种传统优化算法及启发式优化算法。它能够与各种商业软件接口，将专业软件作为各学科或子系统仿真和求解的工具，通过将iSIGHT与商业软件耦合实现数据流的传递，iSIGHT根据分析软件计算结果进行优化，为分析软件提供新的初值，直到计算结果满足设计目标和约束为止。iSIGHT同时是一个开放环境，可以集成用户利用各种计算机高级语言开发的分析计算程序或软件。其他软件ANSYS、I-DEAS等也具有多学科分析优化的功能。

更新设计变量 iSIGHT与仿真分析软件耦合协议 设计任务 分析与分解 输入设计 变量初值 N 满足设计 要求吗？ Y 输出结果 仿真分析软件 输出计算结果 iSIGHT与仿真分析软件耦合协议 iSIGHT优化分析

iSIGHT集成设计优化工作流程

本课程内容

线性规划

一维搜索方法：进退法，黄金分割法，拉格朗日插值，插值和拟合，一元及多元非线性方程求根

无约束优化问题 约束优化问题

多目标函数优化设计 MATLAB优化实现 本课程要求

1、了解优化算法特点。

2、掌握MATLAB软件的优化工具箱函数。

3、熟练掌握MATLAB软件的优化问题的设计分析。

4、注重培养学生的独立设计及开发能力，采用理论与实践相结合，理论讲述与实例分析相结合的方法进行教学，培养和提高学生分析问题和解决问题的能力。 考核方式：笔试（闭卷）

本课程为考试课。

本课程根据教学内容安排3-5次上机作业。

考核主要由四方面综合评定，即平时出勤占10% ，平时作业占20% ，上机实验综合成绩占50% ，随堂考试成绩占20% 。

理论课程主要采用课件授课方式。

第2章 优化设计的数学基础 2.1 向量与矩阵的范数 1. 向量的范数

向量是既有大小又有方向的量，一般用坐标分量来表示，如x?[x1x2...xn]T。向量

的大小用某种数量来表示，称为范数，其中向量的长度就是常用的一种范数，记作x2，即

22x2?x12?x2?...?xn （2.1）

任意向量x?[x1x2...xn]T的范数记为x，它是一个实数，且满足下列3个条件：

nn)，x?0时，x?0；（1）x?0,?x?R(R是n维实向量线性空间

（2）?x??x,???R,x?Rn; ?x,y?Rn.

（3）x?y?x?y,常用的向量范数有：

向量的无穷范数

x??maxxi（或称为l∞范数）

1?i?n向量的1范数

x1??xi（或称为l1范数）

i?1n向量的2范数

x2?(?xi2)2（或称为l2范数）

i?1n1x=[1 2 -3 -8 -9 5] 9 28 1 4 9 64 81 25 2 矩阵的范数

对于给定的n阶方阵A，将A?maxx?0Ax （2.2） x称为矩阵A的范数。

常用的矩阵范数有下面几种：

（1） 矩阵的无穷范数或称为最大行范数

A??max?aij

1?j?nj?1n（2） 矩阵的1范数或称为最大列范数

A1?max?aij

1?j?ni?1n（3） 矩阵的2范数

A2?[?max(ATA)][1 2 3 4]

12T

其中，?max表示矩阵AA的最大特征值。2.2 方向导数与梯度

1 方向导数

首先回顾偏导数的概念，对二元函数，有

(0)(0)?ff(x1(0)??x1,x2)?f(x1(0),x2)|x(0)?lim

?x1?0?x1?x1(0)(0)f(x1(0),x2??x2)?f(x1(0),x2)?f|x(0)?lim

?x2?0?x2?x2根据图2.1，二元函数的方向导数定义为

(0)(0)f(x1(0)??x1,x2??x2)?f(x1(0),x2)?f|x(0)?lim??0?d?(0)(0)f(x1(0)??x1,x2)?f(x1(0),x2)?x1?lim???0?x1??lim?f(x(0)1??x1,x(0)2??0??x2)?f(x?x2(0)1??x1,x)?x2?(0)2 （2.3）

??f?f|x(0)cos?1?|(0)cos?2?x1?x2x1

式中，cosα和cosα

2

分别为d方向与坐标轴x1，x2方向之间夹角的余弦；

??(?x)2?(?y)2。

x2 x(1) d x(0) ρ α2 Δx2 α1 Δx1 0 x1 图2.1 函数的方向导数

n元函数f（x1,x2,…,xn）在点x(0)=[x1(0) x2(0) … xn(0)]T处沿d方向的方向导数为

?f|x(0)??dn?f?f?f|x(0)cos?1?|x(0)cos?2?...?|x(0)cos?n?x1?x2?xn???f|x(0)cos?i?xi?1i （2.4）

2 梯度

（1）二元函数的梯度 根据式（2.3），称向量

?f(x(0))?[?f?x1?fT](0)?x2x??f???x???1? ??f????x2??x(0)为二元函数f（x1,x2）在点x(0)处的梯度。

（2）多元函数的梯度

多元函数f（x1,x2,…,xn）在点x(0)处的梯度为

?f(x(0))?[?f?x1?f?x2...?fT](0) ?xnx?cos?1??cos??2?设d??

?????cos?n??为d方向的单位向量，则沿d的方向导数为

?f|x(0)??f(x(0))Td??f(x(0))cos(?f,d) ?d其中，?f(x(0))代表梯度向量?f(x(0))的范数；cos(?f,d)表示梯度向量与d方向夹角的

余弦。在x(0)处函数沿各方向的方向导数是不同的，它随变化，即随cos(?f,d)所取方向的不同而变化。其最大值发生在cos(?f,d)取值为1时，也就是当d方向和梯度方向重合时其值最大。可见，梯度方向是函数值变化最快的方向，而梯度的范数就是函数变化率的最大值。

例2.1 求目标函数f(x)=x12+x1x2 x(0)=[1 2]T处的梯度。 d a1=30°, a2=60°方向导数 解：

?f?2x1?x2|(1,2)?4,?x1f(x)在x(0)点处的梯度为

?f?x1|(1,2)?1 ?x2??f???x??4?(0)?f(x)??1????

??f??1????x2??x(0)f(x)在x(0)点处梯度的模为

?f(x(0))?42?12?17

2.3 函数的泰勒级数展开

一元函数f(x)的泰勒级数展开，设f(x)在开区间（a，b）内具有直到（n+1）阶的导数，x0

是区间（a，b）内一点，x是以x0为中心的某个领域内的点，则f(x)在x0处的泰勒级数展开为

f(x)?f(x0)?f?(x0)f??(x0)(x?x0)?(x?x0)2?0(?x2) 1!2!1f??(x0)(?x)2 2取级数的前3项近似目标函数，则f(x)可近似表示为

f(x)?f(x0)?f?(x0)?x?对于n元函数f(x)=f(x1,x2,…,xn)，设f(x)在x(0)点的某一领域内存在连续偏导数，则在这个领域内函数f(x)可展开成泰勒级数。即

f(x)?f(x(0))??f?f(0)|x(0)(x1?x1(0))?|x(0)(x2?x2)?...?x1?x2?1??2f?2f?2f(0)2(0)(0)(0)??2|x(0)(x1?x1)?|x(0)(x1?x1)(x2?x2)?|x(0)(x1?x1(0))(x3?x3)???2??x1?x1?x2?x1?x3?1?3f?3f(0)3(0)(0)?[3|x(0)(x1?x1)?|x(0)(x1?x1(0))(x2?x2)(x3?x3)3?x1?x1?x2?x3?3f(0)(0)?|x(0)(x1?x1(0))(x2?x2)(x4?x4)??]???x1?x2?x4

如果忽略二阶以上的各阶小量，则函数可近似表述为

1f(x)?f(x(0))?[?f(x(0))]T(x?x(0))?(x?x(0))TG(x(0))(x?x(0)) （2.5）

2其中，x?[x1x2?xn]T （2.6）

?f(x(0))?[?f?x1?f?x2??fT]x(0) （2.7） ?xn??2f??x2?21??fG(x(0))??2f(x(0))???x2?x1???2??f??x?x?n1?2f?x1?x2?2f2?x2??2f?xn?x2?2f???x1?xn???2f???x2?xn? （2.8）

????2f??2??xn?x(0)

G(x(0))称为f(x)在点x(0)处的Hessian（黑塞）矩阵。 对二元函数，式2.5表示为

f(x1x2)?f(x1(0)(0)x2)??f?f(0)|x(0)(x1?x1(0))?|x(0)(x2?x2)?x1?x21??2f?2f?2f(0)2(0)(0)(0)2???2|x(0)(x1?x1)?2|x(0)(x1?x1)(x2?x2)?2|x(0)(x2?x2)?2??x1?x1?x2?x2?以向量形式表示二元函数f(x)在点x(0)处展开成泰勒二次多项式为

??f?f1??2f?2f?2f2f(x)?f(x)?dx1?dx2??2(dx1)?2dx1dx2?2(dx2)2?

?x1?x22??x1?x1?x2?x2?(0)用矩阵形式表示则为

??ff(x)?f(x(0))????x1(0)?f??dx1?1?dx??2?dx1?x2???2???2f??x2dx2??21??f??x?x?21?2f??dx1??x1?x2??? (2.9) ?2f??dx?2?2??x2?其中，dxi?xi?xi,i?1,2. 根据二元函数的梯度定义，

??f?f(x(0))????x1故式2.9可改写成

?f??,?x2?x(0)T??2f??x22(0)?f(x)??21??f??x?x?21T?2f??x1?x2?? 2?f?2??x2?x(0)f(x)?f(x(0))??f(x(0))x1?x1(0)?函数f(x)在点x(0)处的梯度?f(x重要依据。

(0)????1?x?x?21(0)T1??2f(x(0))x1?x1(0)

??)和Hessian矩阵是计算函数极值以及判断极值点性质的

2.4 无约束优化问题的极值条件

函数的极值点（如极小值点）定义为：对函数定义域内任意一点x(0)，若总存在

f(x(0)??)?f(x(0))，则称x(0)为函数f(x)在x(0)处的局部极小值点。若函数f(x)是单峰函

数，则该极小值点也是f(x)的全局极小值点。因此，求函数的一阶导数（或一阶偏导数）并令其为0是求无约束优化问题函数极值的基本方法。以一元函数为例，如果函数f(x)在某点处一阶导数为0，则该点有可能是函数的极值点，进一步的判断要用到函数f(x)在驻点处的二阶导数f??(x0)。也就是对单调、连续、可微的一元函数f(x)，x=x0为其极值点的充分必要条件为

f?(x0)?0,??0(极小值)?f??(x0)??? （2.10）

0??（极大值）?对于多元函数来说，由式2.7可以看到，如果x(0)是函数f(x)的极小值点，则必有

f(x)>f(x(0))

其中，x是以x(0)为中心的某领域内的点。利用函数取得极值的必要条件，有

[?f(x(0))]?0

因此，可得出下面的不等式：

T1x1?x1(0)??2f(x(0))x1?x1(0)?0 2f(x)?f(x(0))?或

????f(x)?f(x(0))?T1x1?x1(0)G(x(0))x1?x1(0)?0 (2.11) 2????将式2.11应用到一元函数的情形，可知，这也就是式2.11。 对于对称矩阵G，若存在非零向量x，满足

xTGx?0 （2.12）

则称G为正定矩阵。对于极大值问题，则要求G为负定矩阵（即），因此，n元函数f(x)在点x(0)处取得极值的充分必要条件是

??f?f(x(0))????x1?f...?x2?f?T?[00...0]?0 (2.13) ??xn?T??2f??x2?21??f(0)2(0)G(x)??f(x)???x2?x1???2??f??x?x?n1?2f?x1?x2?2f2?x2??2f?xn?x2?2f???x1?xn??2?f???x2?xn? （2.14）

????2f??2??xn?x(0)（若为正定，则为极小值；若为负定，则为极大值。）

对于二元函数y=f(x1,x2), x(0)为其极值点的充要条件可表示为

??f(0)?f(x)????x1?f??0（必要条件）

?x2?(0)?xT??2f??x2(0)G(x)??21??f??x?x?21?2f??x1?x2??正定或负定（充分条件） ?2f?2??x2?x(0)?2f?2f(0)(0)(0)(0)且2?0时，(x1,x2)为极大值点，2?0时，(x1,x2)为极小值点。 ?x1?x1对称矩阵G的正定性除可按定义式2.10判定外，更实用的方法是利用矩阵的各阶主子式的正负号进行判定。若G的各阶主子式均大于0，则G为正定矩阵；若G的各阶主子式的正负号负、正相间，即奇数阶主子式小于0，偶数阶主子式大于0，则G为负定矩阵。

?1021???例2.2 判断矩阵A?253的正定性。 ????134??解：

因为a11>0，则

即矩阵A的各阶主子式的值均大于0，故矩阵A为正定。 例2.3 判断矩阵

??12?2??

A???35?3????0?1?2??的正定性。 clc;

A=[-1 2 -2;-3 5 -3;0 -1 -2]; for i=1:3

A1(i)=det(A(1:i,1:i)); end -1 1 -5

2x12x2例2.4 试求函数f(x)?的极值。 ?22[0 0]T

G=[1 0

0 1] 正定 极小值

2.5 凸集与凸函数

经典优化算法大多属于沿某一搜索方向的局部优化算法，要求目标函数和约束条件满足一定要求，如要求目标函数及约束条件为单峰函数或单调函数，或者说要求目标函数和约束函数为凸函数，对应的解集为凸集。相对于全局优化算法来说，凸集、凸函数的概念对于局部优化算法更为重要。如果已知目标函数为凸函数且求解域为凸集，则求解出的局部极值点也是全局极值点，否则求出的局部极值点就不一定是全局极值点。对于约束优化问题应验证解是否满足约束条件。

1 凸集

设集合S?R，如果对于任意的x1，x2∈Rn，有

n?x1?(1??)x2?S,???[0,1] (2.15)

则称S是一个凸集。该定义用文字描述就是：对于一个点集（或区域），如果连续其中任意两点x1，x2的线段都全部包含在该集合内，就称该点集为凸集，否则为非凸集。 反过来，如果是一个凸集，点x1,x2,…,xm∈S，则这m个点的组合为凸集：

??x?S

iii?1m其中，

??i?1mi?1,?i?0（i=1，2，…，m）

例2.5 证明超平面H?x?R|px?c是一个凸集。其中p∈Rn是超平面的法向量，且不为0；c是一个标量。

2 凸函数

定义：设f(x)为定义在非空凸集x?R上的实值函数，若对任意x1,x2∈Rn以及数??[0,1]，恒有

n?nT?f[?x1?(1??)x2]?(?)?f(x1)?(1??)f(x2)

则称f(x)为x上的下凸（上凸）函数。 若对任意x1,x2∈Rn以及数??[0,1]，恒有

f[?x1?(1??)x2]?(?)?f(x1)?(1??)f(x2)

则称f(x)为x上的严格下凸（上凸）函数。 性质：

（1）设f(x)为定义在凸集上Rn上的凸函数，则对于任意实数β≥0，函数βf(x)也是定义在Rn上的凸函数。 （2）设f1(x)和f2(x)为定义在凸集Rn上的两个凸函数，α，β为不同时为0的实数，则f(x)= αf1(x)+βf2(x)仍然是定义在Rn上的凸函数。

（3）设f(x)为定义在凸集Rn上的凸函数，则对于任意实数β，集合S?x|x?R,f(x)???n?是凸集。

（4）设f(x)为定义在凸集Rn上的可微凸函数，若存在点x*∈Rn，使得对于所有的x∈Rn有[?f(x)][x?x]?0，则x*是在f(x)上的最小点（全局极小点）。该性质说明，函数f(x)在极值点x*处的梯度与曲线上两点割线构成的矢量的夹角α≤π/2，在极值点处函数的梯度与过极值点的切线垂直。 凸规划

*T*

?a110??x1??b1???? ???a??211??x5??b2?选a11为消元主元，结果为

b1????10??x1??a11? ?01??x???b1???5??b2?a21??a11???从而得

x1?b1,a11x5?b2?a21b1 (3.4) a11用x1替换x5的约束方程为

?1a11??x4??b1??0a??x???b? (3.5)

21??1???2?选a21为消元主元，结果为

b2??b?a?10??x4??111a21?? (3.6) ?01??x???b2????1????a21??从而得

x4?b1?a11b2,a21x1?b2 (3.7) a21观察式（3.4）和式（3.7），注意到当

?b??5x1?min?i??min??1?ai1?3?b23??3 (3.8) ??1?a211时，可得到基本可行解，否则x1的取值不能使所有基本变量的取值成为基本可行解，如式3.4。出基变量与式3.8有着必然联系。式3.8中最小比值对应的右端项元素或系数行下标对应的基本变量就是出基变量，且比值分母系数下标就是进基变量对应的主消元素。验证式3.3可知，变量x5为被替换出的变量，与计算机结果一致。由此可见，目标函数与约束条件配合，约束方程求解完成后，将基本变量结果代入目标函数中，通过比较非基本变量系数大小决定进基变量（此时认为非基本变量取值不为0）。对求最大值问题，系数较大的非基本变量为进基变量；对求最小值问题，系数较小的非基本变量为进基变量。确定了进基变量后，则计算右端项与进基变量在约束方程中列系数向量对应元素的比值，比值小者的行下标对应的基本变量为出基变量。

将式3.1中的约束矩阵A按列表示： A=[ p1 p2 … pn]

pi (i=1,2, …,n)是m×1的向量，pi的分量是各约束方程中与设计变量xi对应的系数。设矩阵A的秩为m将设计变量x分解为基本变量和非基本变量，即x=[ xE xN ]T，相应地，系数矩

阵A也分解为两部分，一部分为基矩阵，用E表示；另一部分为非基矩阵，用N表示。于是A=[ E N ]，E=[ p1 p2 … pm]，N=[ pm+1 pm+2 … pn]。式3.1中的约束条件重新表示为

?E即

?x?N??E??b （3.9） ?xN?ExE?NxN?b

上式两端左乘E-1并移项，得

xE?E?1b?E?1NxN （3.10）

xN的分量为自由变量，它们取不同的值，就会得到方程组的不同解。特别地，如果令xN=0，则得到

?xE??E?1b?x?????? （3.11）

?xN??0?称解x为方程组Ax=b的一个基本解。如果E-1b≥0，则称

?xE??E?1b?x??????

x?N??0?为满足约束条件的基本可行解（非负的基本解）。 2 基本可行解的转换

基本可行解的转换就是得到某组基本可行解后，如何进一步改进获得更好的解，最终达到最优基本可行解。从上节的分析可知，获得改进的解是在满足约束条件的前提下使目标函数值增加（对于最大值问题）或使目标函数值减少（对于最小值问题）的解。设已经获得一组基本可行解：

x(0)?E?1b????

0??其对应的目标函数值记为f0。设x???xE??是另一组基本可行解，将其代入目标函数中，得 x?N?Tf?cTx?cE?xE?cTN???xN?TT?1?1TT?1TT?1 ?cExE?cTNxN?cE(Eb?ENxN)?cNxN?cEEb?(cN?cEEN)xN （3.12）

???f0?j?m?1?(cnj?cEpj)xj?f0?TE?1j?m?1?(cnj?zj)xj如果每次将一个非基本变量转换为基本变量，并假定目标函数的最大值，那么由式3.12可知，选择求和项中系数cj?zj>0最大者对应的非基本变量进入基本变量，当该变量取正值时，可使目标函数增加最多。如果所有非基本变量的系数cj?zj不再有正值，则说明非基

本变量的进基运算不会再使目标函数值增加，此时就终止计算，输出最优结果。非基本变量转换为基本变量后，则要将上一轮迭代中的一个基本变量转换为非基本变量，判断出基本变量表的条件由约束方程的消元计算格式给出。

设xk为进基变量，不失一般性，认为xk替换xr，在约束方程中相应地用pk代替pr，

?kpk??a1?k?ank??T,pr??a1?ra2?r?anr??,新的基矩阵为E（1）a2。为表示清楚，

T(k)pk将标记为pr(k)，其元素相应地表示为pr表示为

?k?a1?(k)?ka2(k)??ank(k)T?。因此基矩阵E

（1）

E(1)?p1?p2?pr(k)?pm??r(k)?10?a1??(rk)?01?a2???????(k)????arr??????(k)??00?amr?1?x?E?E(1)b (3.13)

?0???0??? ??0?????1??变量xk由0变为正值后，新方程组的解为

其中，基矩阵E（1）的逆矩阵为

??r(k)a1?10??a?(k)rr??(rk)?01??a2??(k)arr?????1E(1)???????1??(k)arr??????(k)?00??amr??(k)arr???0???0????? (3.14) ?0??????1???初始基矩阵为单位矩阵，在基矩阵中用pk将标记为pr，就是将pk放在第r列的位置上。对

于式3.14所示矩阵，其逆矩阵只有第r列的元素发生变化，对角线上的元素为相应元素的倒数，该列上其他元素等于原来的元素取反再除以对角线上的元素。式3.13的解为

xr(k)?br?（3.15） ?(k)arrx(k)i?(k)br?air?(k)xr(k),i?1,2,?,m,i?r（3.16） ?bi??(k)?bi??air?arr由于设计变量非负，并假设设计变量在约束方程中的系数和右端项元素均大于0，因此xr(k)

≥0自然满足；而式3.16中的变量非负意味着

?(k)br?airbi??(k)?0

?arr即

bi?br???0（3.17） (k)(k)??airarr当

??xk??bi??br?(k)??min,a?0?(k)ir?（3.18） (k)??arr?air?时，则能保证用xk替换xr，使所有基本变量非负。在线性规划中，通常称

?j?cj?zj（3.19）

j为判别数或检验数。对极小化问题，进基变量xk的下标与?k?ck?zk?min(cj?zj)项的

下标k对应，非基本变量xk的引入使目标函数f?f0?j?m?1??njxj增加最多（最大值问题），

或减少最多（最小值问题）。出基变量xEr的下标与???b??br???0?项的下标?min?i,aik“r”

??arka?ik?对应。式3.18和式3.19是单纯形法求解线性规划问题的重要判别式和检验规则。

3 单纯形法的计算步骤

例3.2 用单纯形法求下述问题的最优解：

minz??6x1?3x2?2x3s.t.2x1?x2?x3?6x1?x3?2x1,x2,x3?0解：

首先引入松弛变量x4，x5，把数学模型化为标准形式：

minz??6x1?3x2?2x3?0x4?0x5s.t.2x1?x2?x3?x4?0x5?6x1?0x2?x3?0x4?x5?2x1,x2,x3,x4,x5?0引入松弛变量后，变量总数n=5，约束方程数m=2（不含变量大于等于0的约束），因此有n-m=3个非基本变量，并令其为0进行求解。一个简单的做法就是选松弛变量x4，x5为基本变量，因此其系数列向量为单位基向量E??p4和非基本变量表示，得 x4=6-2 x1- x2- x3

?10?p5????。将基本向量用在端列向量

01??x5=2- x1 - x3

令非基本变量x1=x2= x3=0，则基本解为 x4=6，x5=2

记初始基本解为 x(0)=[ 0 0 0 6 2 ]T

此解又满足式，xj≥0，j=1，2，…，5，因此，x(0)是基本可行解，它对应着凸可行域的一个顶点。一般来说，对于约束条件全为“≤bi”形式（i=1，2，…，m）的线性规划问题，通过引入松弛变量，可找到一个初始基本可行解。

将初始基本可行解x(0)=[ 0 0 0 6 2 ]T代入目标函数表达式中，得 z=-6×0+3×0-2×0+0×6+0×2=0

非基本变量x1=x2= x3=0对应的目标函数值z=0，为找出进基变量，将x4，x5代入原式得

z??6x1?3x2?2x3?c4(6?2x1?x2?x3)?c5(2?x1?x3)?6c4?2c5?(?6?2?c4?c5)x1?(3?c4)x2?(?2?c4?c5)x3?x1??10211??????

??c4c5????101???x2?01??????x??3?T?1?cTEE?1b?(cTN?cEEN)xN??6x1?3x2?2x3?10??6??c5????2?????63?2???c401?????由上式看出，x1增加一个单元，目标函数值z下降6个单位；x2增加一个单元，目标函数值

z增加3个单位；x3增加一个单元，目标函数值z下降2个单位。如果使x1和x3成为正值，都能使目标函数向极小改善，故当前解不是最优解。

选择的进基变量应是使目标函数改善最大的非基本变量进基。在目标函数中，应选择有负系数，且负系数绝对值最大的非基本变量进基。例3.2中，目标函数中x1的系数为-6，，所以x1应选为进基变量。这是因为当x1由当前的0变为正值时，使目标函数下降最大。

出基变量的选择也不是任意的，原则是要保证变量满足非负条件。这两个方程均有变量x1要从0上升为正值，要从x4和x5这两个变量中，选择一个出基（从当前值下降到0），而另两个非基本变量x2和x3还要保持为0. 由式

x4?6?2x1?x2?x3x5?2?x1?x3得

x4?6?2x1x5?2?x1

x1从0变为正值时，x4和x5可从当前值下降为0，但不能成为负值（因为变量要满足非负条件）。由上式可以看出，x1取x1??min?min(bi2)??2时，x4和x5均为非负（x4=2>0，ai11x5=0），且θmin使x5=0，因此选x5为出基变量。

已经确定了x1进基变量和出基变量x5，就可进行新的消元求解。继续求解：

2x1?x2?x3?x4?0x5?6x1?0x2?x3?0x4?x5?2

此时，令非基本变量，x2= x3= x5=0，求x1和x4的解。用初等变换（消元）求解上面的方程，结果为

-2 x5+ x2- x3+ x4+ 0x1=2 x5+0 x2+ x3+0 x4+ x1=2 解得

x1=2，x4=2

再将x1和x4用约束方程的非基本变量表示，得 x4=2- x2+ x3+2 x5 x1=2- x3- x5

将其代入目标函数，得

z=-6（2- x3 – x5）+3 x2-2 x3+c4（2- x2+ x3+ 2x5）+ c5 x5 整理得

z=-6×2+ c4×2+（3- c4）x2+（-2+6+ c4）x3+（c5+6+2 c4）x5 由于所有非基本变量的系数均为非负，没有进一步可使目标函数减少的非基本变量，因此迭代结束，最优解为x1=2，x2=0，x3=0，f=-12。 单纯形法的一般步骤可归纳如下：

（1）把线性规划问题的约束方程组表达成标准形方程组，然后找出一个基本可行解作为初始基本可行解。对等式约束或“≥bi”形式的约束，如果不易得出初始基本可行解，则需用辅助方法得出，如大M法、两段法等。

（2）若基本可行解不存在，即约束条件有矛盾，则问题无解。

（3）若基本可行解存在，则从初始基本可行解作为起点，根据θ规则（式3.18）和判别数ζ（式3.19），引入非基本变量取代某一基本变量，找到使目标函数值更优的另一基本可行解。 （4）按步骤（3）进行迭代，直到对应检验数满足最优性条件（这时目标函数值不能再改善），即得到问题的最优解。

（5）若迭代过程中发现问题的目标函数值无界，则终止迭代。

4 单纯形法列表计算

用单纯形法求解线性规划问题的过程，实际上就是解线性方程组。只是在每次迭代中，要按一定规则选择自由未知量，以便能够得到改进的基本可行解。这个求解过程可以通过变换单纯形法来实现。单纯形是以约束方程的增广矩阵为中心构造的一种变换表格，见表3.1。 表3.1 单纯形表 设计变量 基本变量 xn-m+1 xn-m+2 . . . xn ζj=cj-zj 非基本变量 cf ch 0 0 . . . 0 x1 c1 p1 x2 c2 p2 … … … xn-m cn-m pn-m xn-m+1 0 pn-m+1 xn-m+2 0 pn-m+2 … … … xn 0 pn b1 b2 . . . bm f(x) b bi/ air 应用单纯形表求解线性规划问题的基本步骤如下：

（1） 将线性规划问题写成标准形式，在单纯形表中对应位置填入目标函数变量的系数、

约束矩阵中变量的系数和右端项，将增加的m个松弛变量作为基本变量，初始化单

纯形表；

（2） 计算minζj=cj-zj，确定进基变量xL，计算minθi= bi/ bir，确定出基变量xi；

（3） 以L列、i行对应的约束矩阵中的元素aiL为主元进行消元，对矩阵元素和右端项更

新，并计算目标函数值；

（4） 重复第（2）和（3）步，直到最后一行各列ζj≥0时消元中止，此时与基本变量对

应的右端项bi就是所求的最优解。

例3.3 用单纯形法求解下面的线性规划问题：

f(x)??2x1?x2?x3s.t.3x1?x2?x3?60x1?x2?2x3?10 x1?x2?x3?20x1,x2,x3?0解：

线性规划的标准形式为

minf(x)??2x1?x2?x3?0x4?0x5?0x6s.t.3x1?x2?x3?x4?0x5?0x6?60x1?x2?2x3?0x4?x5?0x6?10 x1?x2?x3?0x4?0x5?x6?20x1,x2,x3,x4,x5,x6?0（1） 构造初始单纯形表 初始单纯形表如表3.2。 表3.2 初始单纯形表

设计变量 基本变量 x4 x5 x6 ζj=cj-zj 非基本变量 cf Cb 0 0 0 x1 -2 p1 3 1 1 -2 x2 -1 p2 1 -1 1 -1 x3 1 p3 1 2 -1 1 x4 0 p4 1 0 0 0 x5 0 p5 0 1 0 0 x6 0 p6 0 0 1 0 60 10 20 f=0 60/3 10/1 20/1 b bi/ air

（2） 消元

此问题是求目标函数的最小值，判别式minζj=cj-zj为求各列的最小值，L=1，再计算minθi= bi/ bir，得i=2。因此x1进基，x5出基。得出新的单纯形表，并以a21为主元进行消元，结果见表3.3。

表3.3第1次消元结果 设计变量 基本变量 x4 x1 x6 ζj=cj-zj 非基本变量 cf ch 0 -2 0 x1 -2 p1 0 1 0 0 x2 -1 p2 4 -1 2 -3 x3 1 p3 -5 2 -3 5 x4 1 p4 1 0 0 0 x5 0 p5 -3 1 -1 2 x6 0 p6 0 0 1 0 30 10 10 f=-20 30/4 10/-1 10/2 b bi/ air

由表3.3可知，minζj=cj-zj对应的最小列号为L=2，相应地minθi= bi/ air对应的行号为i=3，因此为x2进基变量，x6为出基变量，得到新的单纯形表，并以a32为主元进行消元，结果见表3.4。

表3.4 第2次消元结果 设计变量 基本变量 x4 x1 x2 ζj=cj-zj 非基本变量 cf ch 0 -2 -1 x1 -2 p1 0 1 0 0 x2 -1 p2 0 0 1 0 x3 1 p3 1 1/2 -3/2 1/2 x4 1 p4 1 0 0 0 x5 0 p5 -1 1/2 -1/2 1/2 x6 0 p6 -2 1/2 1/2 3/2 10 15 5 f=-35 b bi/ air 经过两次消元，判别式ζj=cj-zj≥0，求解结束，最优解为x*=[ 15 5 0 ]T,f(x*)=-35。

3.4 改进的单纯形法

单纯形法进行消元计算时对约束矩阵的所有元素都进行变换运算，计算量很大，为了减少计算量，提出改进的单纯形法。 1 改进的单纯形法的基本思想 对于给定的线性规划问题

maxf?cTxs.t.Ax?b x?0应用单纯形法求解时，基本变量，目标函数和判别数用下面表达式计算：

xE?E?1b?E?1NxNT?1T?1TT?1f?cEEb?(cTN?cEEN)xN?f0?(cN?cEEN)xN (3.20) TT?1?N?cTN?cEEN由式3.20看出，每个计算式都含有关于E-1的运算，因此E-1是计算基本变量、判别数的关键。对于求最大值问题，如果ζN的分量中有正数，选择一个非基本变量将其转换为基本比变量，相应地用新的可行基的逆矩阵再计算更新的基本变量。目标函数和判别数。 设用非基本变量xk替换基本变量xr，则替换前后的基矩阵为

E??p1~E??p1p2?pr?1p2?pr?1prpkpr?1?pm? pr?1?pm?

E?1E?E?1?p1?E?1p1E?1p2?E?1pr?1E?1pr??e1e2?er?1er?p2?pr?1prpr?1?pm?er?1?em?E?1pr?1?E?1pm?I?

???e~E?1E?E?1?p1p2?pr?1E?1pkpkpr?1?pm??E?1p11E?1p2?E?1pr?1E?1pkE?1pr?1?pm

e2?er?1?1er?1?em??~?1~?1?记Erk?EE，则Erk?EE就是式3.14。 ??E?1E因为Erk???~?1?~~?E?1。如果初始矩阵E为单位矩阵，则?E?1E，因此，E?1?Erk~~?1?E?1?ErkE?1?Erk。对于第l次迭代，基矩阵的逆矩阵为 ~1~?1~?1~?1~?1~?1E(??EE?E?El)(l),rk(l?1),rk(1),rk(l),rkE(l?1) (3.21)

式3.21中，下标rk随每次迭代进基变量和出基变量的下标而改变。用式3.21进行迭代计算时，约束方程中变量的系数、右端项保持不变，每次计算只需保存前一次基矩阵的逆矩阵。 在一次换基迭代过程中，在确定了进基向量与出基向量之后，真正起作用的是进基列向量pk与前一个基矩阵的逆矩阵E-1。

2 改进的单纯形法的计算步骤

设E(0)为线性规划问题的一个初始可行基。步骤为：

(1) 计算E-1并代入式3.19，计算基变量、目标函数和检验数。

TT?1T?1?cTxE?E?1b,fE?cEEb,?NN?cEEN

(2) 对检验数进行判定。若?N?cN?cEEN的分量全为非正（求最大值问题），则计算结

?1?E?1b?束，x???即为最优解，

0??TfE?cEE?1b为最优值。否则转入第（3）步。

(3) 确定进基列向量。令?k?max?j,j?1,2,...,n，则pk为进基列向量。并计算

????E?1pk??a1?kpk?ka2??amk??T，如果所有的ark?≤0(i=1,2,…,m)，则此问题?ark?>0，则转入第（4）步。 无最优解；若至少有一个，如ark??b?br?为主元，确定pr?出基。?min?r,i?1,2,...,m?，以ark(4) 确定出基列向量。令?r?

??ark?ark?(5) 计算E~?1?1?1?取代单元方阵中er得来。 ?ErkE。其中，Erk由pk~?1?1rk?1?1rkE~~?1(6) 计算xE?Eb?EEb?Ex,?N?cN?cEEN，返回第（1）步。 例3.5 用改进的单纯形法计算下面的线性规划问题：

maxf(x)?25x1?90x2,x?R2s.t.2x1?x2?150x1?3x2?270x2?80x1,x23?0解：

将所给问题化成线性规划问题的标准形式：

maxf(x)?25x1?90x2?0x3?0x4?0x5s.t.2x1?x2?x3?150x1?3x2?x4?270x2?x5?80xj?0,j?1,2,...,5?21100?T??设A?13010,c??2590000?,b??150????01001??(1) 构造初始基矩阵 初始基矩阵为E0??p3

27080?

Tp4?100??100??，同时得E?1??010?，所以有

p5???0100???????001???001??xE0?x3??100??150??150???E?1b??010??270???270? ??x40???????????x5???001????80????80???100????000? T?1cE00??0100E0??0????001??(2) 计算各非基本变量的判别数

?2??1??25 T?1???1?c1?cEEp?25?000001????0???1???90 T?1?2?c2?cE00??30E0p2?90??0????1??(3) 确定进基变量

?1?25,?2?90???2，对应非基本变量，确定为x2进基变量。同时计算 根据max?

?100??1??1???1???1?

E0?1p2??010????????001????1????1??(4) 确定出基变量

?bi??150,i?1,2,3??min?根据θ规则，求出?r?min??1?ai2?为出基变量。于是得到新的基矩阵E1??p3

(5) 计算新的基矩阵的逆阵

270180?80|r?3,，x5被确定??1?1p4p2?，相应的cE1??0090?。

?101??10?1????01?3? ?1E1?1?(E0E1)?1??013???????001???001???x3??10?1??150??70???E?1x??01?3??270???30?

xE1??x41E0???????????x2???001????80????80???10?1????0090? T?1cE090??01?31E1??0????001??用E1代替E0，重复步骤(2)~(5)。 (6) 计算各非基本变量的判别数

?1?1?1?10?1??2???1??25 T?1?1?c1?cE090??01?30E1p1?25??0??????001????0???10?1??0??01?1??0???90 T?1???5?c5?cEEp?0?0090115??????001????1??(7) 再确定进基变量

?1?25,?5??90???1，对应非基本变量x1，确定x1为进基变量。同时计算 根据max??10?1??2??2???1???1?

E1?1p1??01?3????????001????0????0??(8) 再确定出基变量

根据θ规则，求出?r?min??bi??70,i?1,2,3??min??2?ai1?30180?30|r?2，x4被确定为??0?1出基变量。于是得到新的基矩阵E2??p3

(9) 计算新的基矩阵的逆阵

p1T2590?。 p2?，相应的cE2??0?121??1?25????01?3? ?1E2??013??????1??001???00?xE2?x1??1?25??150??10???E?1x??01?3??270???30? ??x2E0?2?????????1??x3???00???80????80???1?1?25????02515? T?1cE2590??01?32E2??0???1??00??1?25??0??01?3??1???25 T?1???3?c4?cEEp?0?02590224?????1??00???0???1?25??0???0???15 T?1?5?c5?cE2590??01?32E2p5?0??0?????1??00???1??因判别数ζj均小于0，迭代结束。最优解为x??30801000?,fmax?7950。

T3.5 优化工具箱

用MATLAB优化工具箱求线性规划的函数linprog（）来求解此问题。

线性规划问题数学模型 min CTX s.t. AX≤B AeqX=Beq Lb≤X≤Ub

其中：C，X，B，Beq，Lb，Ub为向量；A，Aeq为矩阵 函数输入变量

X= linprog（C,A,B）[]

X= linprog（C,A,B,Aeq,Beq） [] X= linprog（C,A,B, Aeq,Beq,Lb,Ub） X= linprog（C,A,B ,Aeq,Beq,Lb,Ub,x0）

X= linprog（C,A,B, ,Aeq,Beq,Lb,Ub,x0,options）

函数左端输出变量 [x]=linprog()

[x,fval]= linprog()

[x,fval,exitflag]= linprog()

[x,fval,exitflag,output]= linprog()

[x,fval,exitflag,output,lambda]= linprog()

X返回目标函数最优解,fval返回目标函数最优值,exitflag终止迭代的标志（整数值）,output输出迭代次数,lambda解x的拉格朗日乘子 exitflag=1 exitflag=0 exitflag=-2 exitflag=-3

output, algorithm output, funcCount output,iterations output,message

第4章 一维搜索方法

是指求解一维目标函数f(x)的极小点和极小值的数值迭代方法，可归纳为单变量函数的极小化问题。虽然优化设计中大部分问题是多维问题，一维问题的情况较少，但是一维优化方法是优化方法中最基本的方法，在数值迭代过程中都要进行一维搜索。另外，很多多维优化问题最终归结为一维优化问题来处理。

如果确定了迭代点x(k)及其搜索方向d(k)，那么迭代所得的新点x(k+1)将取决于步长a(k)，即

x(k+1)= x(k)+ a(k) d(k)，k=0，1，2，… (4.1)

由式4.1可知，不同的步长a(k)会得到不同的迭代点和不同的目标函数值f(x(k+1))。一维优化问题的目的是在既定的迭代点x(k)和搜索方向d(k)下寻求最优步长a(k)，使迭代产生的新点x(k+1)的函数值最小，即

min f[x(k)+ a(k) d(k)]

在初始迭代点x(k)和搜索方向d(k)确定后，就把求解多维优化问题的极小值变成求解一个自变量即最优步长a的最优值的一维问题了。即求一元函数 f(x(k+1))=f(x(k)+ a d(k))=φ( a) (4.2) 的极值问题。

一维搜索的优化方法很多，常用进退法、黄金分割法和二次插值法。 一维搜索方法的步骤：

(1) 确定初始搜索区间[a,b]，即最优步长a所在的区间[a,b]。搜索区间应为单峰区间，并且在区间内目标函数应只有一个极小值。

(2) 在搜索区间[a,b]内寻找最优步长a，使目标函数式4.2达到极小值。

4.1 确定初始单峰区间的方法—进退法

原理

进退法也称为外推法，是一种通过比较函数值大小来确定单峰区间的方法。由单峰函数的性质可知，在极小点左边函数值应严格下降，而在极小点右边函数值应严格上升。因此，从某一种给定的初始点x0出发，以初始步长a沿着目标函数值的下降方向，逐步前进（或后退），直至找到相继的3个计算点的函数值出现“大-小-大”的趋势为止。 利用进退法确定搜索区间[a,b]的步骤如下： (1) 任取x0，步长a>0，取x1=x0+a。 (2) 后退运算。φ(x1)>φ(x0)，则令a=2a（步长加倍），x2=x0-a。分以下两种情况：若φ(x2)<φ(x0)，则令x1=x0，x0=x2。重复(2)；若φ(x2)>φ(x0)，则停止a=x2，b=x1。 (3) 前进运算，若φ(x1)<φ(x0)，则令a=2a（步长加倍），x2=x1+a。若φ(x2)<φ(x1)，则令x0=x1，x1=x2，重复(3)；若φ(x2)>φ(x1)，则停止，a=x0，b=x2。 程序框图

4.2 黄金分割法

1 基本原理

又称为0.618法，它通过不断缩短搜索区间的长度来寻求一维函数f(x)的极小点。对于单峰函数f(x)，在其极值存在的某个区间[a,b]内取若干点，计算这些点的函数值并进行比较，总可以找到极值存在的更小区间。在这更小区间内增加计算点，又可以讲区间进一步缩小。当

区间足够小，即满足精度要求时，就可以用该区间内任意一点的函数值来近似表达函数的极值。

设单变量函数f(x)在区间[a,b]上有定义，若存在一点x*(a b)，使得f(x)在区间[a,x*]上严格单调减，f(x)在区间[x*,b]上严格单调增，则称f(x)是区间[a,b]上的（下）单峰函数。显然x*是f(x)在区间[a,b]上的唯一的极小值点。

根据（下）单峰函数所具有的性质，对在某区间[a,b]上的（下）单峰函数f(x)可采用黄金分割法搜索其在区间[a,b]内得极小值点。 2 计算方法

设区间[a,b]的长度为L，在区间内取点λ1，将区间分割为两部分，线段aλ1的长度记作λ，并满足

?L?L????q2且2??L??

?1?5，取正根2由上式有??L??L?0，两边同除L2，得q2+q-1=0，则有q?q?5?1?0.6180339887。q称为区间收缩率，它表示每次缩小所得的新区间长度与缩2小前区间长度之比。

在一维搜索时，在区间内取两对称点λ1和λ2，并满足

q??2L??1L??2??0.618 ?2?2显然，经一次分割后，所保留的极值存在的区间要么是[a, λ2]，要么是[λ1，b]。而经k次分割后，所保留的区间的长度为??qL?(0.618)L。

由于区间收缩率q是一个近似值，每次分割必定带来一定的舍入误差，因此，分割次数太多

时计算会失真。经验表明，黄金分割的次数k应限制在11以内。

kkk4.3 拉格朗日插值多项式

是一种显式公式,它将pn(x)表示为一组插值基函数的线性组合。 1 线性插值

设函数y=f(x)在给定的互异节点x0，x1上的函数值分别为y0=f(x0)，y1=f(x1)，若能构造一个函数

p1(x)=a+bx (4.3)

使它满足p1(x0)=y0，p1(x1)=y1，则式4.3所示的插值问题称为线性插值。

线性插值的几何意义是过曲线y=f(x)上的两点(x0,y0)和(x1,y1)作一直线，用p1(x)近似值代f(x)。

对于给定的两点(x0,y0)和(x1,y1)，某一点x(x0

p1(x)?y0?记

y1?y0x?x1x?x0(x?x0)?y0?y1

x1?x0x0?x1x1?x0

l0(x)?x?x1x?x0 ,l1(x)?x0?x1x1?x0l0 (x)和 l1 (x)称为线性插值基函数。

l0 (x)和 l1 (x)实质上是满足条件

l0 (x0)=1, l0 (x1)=0; l1 (x0)=0, l1 (x1)=1的一次插值多项式。

线性插值的解可以表示为插值基函数l0 (x)和 l1 (x)的线性组合，其组合系数为y0和y1,即 p1(x)= y0 l0 (x)+ y1 l1 (x) 2 二次函数插值

又称为抛物线法，基本思路为：利用目标函数在若干点的信息和函数值，构造一个与目标函数相接近的低次插值多项式，然后求该多项式的最优解作为原函数的近似最优解。随着区间的逐次缩小，多项式最优点与原函数最优点之间的距离逐渐缩小，直到满足一定精度要求时终止迭代。 基本原理

设目标函数f(x)在点x1,x2,x3(x1< x2< x3)上的函数值分别为f1=f(x1),f2=(x2),f3=(x3),且满足f1> f2和f2< f3即满足函数值呈“大-小-大”的趋势。于是可通过原函数曲线上的3个点p1(x1,f1), p2(x2,f2), p3(x3,f3)做一条二次曲线，此二次插值多项式为p(x)=a+bx+cx2。 为了求解p(x)的极小值，对p(x)求导数，并令其为0，即

p?(x)?b?2cx?0

解得二次函数极小点

xp??b (4.5) 2c为求得，应求出式4.5中的待定参数b和c。

根据插值条件，插值函数p(x)与原函数f(x)在插值点p1, p2, p3处函数值相等，得 p(x1)=a+bx1+cx12=f1 p(x2)=a+bx2+cx22=f2 p(x3)=a+bx3+cx32=f3 求得

a?(x3?x2)x2x3f1?(x1?x3)x1x3f2?(x2?x1)x1x2f3(x1?x2)(x2?x3)(x3?x1) (4.6)

2222(x2?x3)f1?(x3?x12)f2?(x12?x2)f3b?(x1?x2)(x2?x3)(x3?x1)c?(x2?x3)f1?(x3?x1)f2?(x1?x2)f3(x1?x2)(x2?x3)(x3?x1)把式4.6代入式4.5，既得二次插值函数极小点的计算公式：

22221?(x2?x3)f1?(x3?x12)f2?(x12?x2)f3?xp??? (4.7)

2?(x2?x3)f1?(x3?x1)f2?(x1?x2)f3?为便于计算，将式4.7改写为

xp?0.5(x1?x3?C1) (4.8) C2式中，

C1?f3?f1(f?f)/(x2?x1)?C1 ,C2?21x3?x1x2?x3将只用一次二次插值计算所得xp的作为函数的极小点x*的近似解往往达不到精度要求，为此需缩短区间，进行多次插值计算，使xp不断逼近原函数的极小点x*。 迭代过程

根据原区间中x2和xp的相对位置f2和fp函数值和的大小，区间缩短有4种情况： 若xp> x2，f2>fp，则以[x2, x3]为新区间，令x1←x2, x2←xp, x3不变； 若xp> x2，f2≤fp，则以[x1, xp]为新区间，令x3←xp, x1, x2不变； 若xp≤x2，f2> fp，则以[x1, x2]为新区间，令x3←x2, x2←xp, x1不变； 若xp≤x2，f2≤fp，则以[xp, x3]为新区间，令x1←xp, x2和 x3不变。 3 n次拉格朗日插值多项式 插值多项式的存在唯一性

所谓构造插值多项式就是要构造一个不超过n次的多项式

pn(x)?a0?a1x???anxn

使其满足以下插值条件：

a0?a1x0???anx0?y0a0?a1x1???anx1?y1?a0?a1xn???anxn?yn式4.9是一个关于待定参数a0，a1，… ，an的n+1阶线性方程组，其系数矩阵行列式

n1x0?x0nnn (4.9)

Vn?1x1?x1n???n1xn?xn?0?j?i?n?(x?x)

ij称为Vandermonde(范德蒙)行列式。当时xi≠xj(i≠j)，满足条件(4.9)的插值多项式存在且唯一。

n次拉格朗日插值多项式

直接通过解线性方程组(4.9)来确定插值多项式的系数，一般计算工作量较大，且得出的多项式不便于应用。因此常用其他一些方法来构造插值多项式。下面采用所谓插值基函数的方法来构造一种具有一定特点的插值多项式。

对于给定的n+1个互异点xi (i=0,1,…,n)处的数据f(xi)=yi (i=0,1,…,n)，令

nx?xj(x?x0)(x?x1)?(x?xi?1)(x?xi?1)?(x?xn)li(x)???,i?0,1,...,n(4.10)

(xi?x0)(xi?x1)?(xi?xi?1)(xi?xi?1)?(xi?xn)j?0xi?xjj?i称li(x)为关于节点xi的n次插值基函数，且li(x)满足

?0,j?i?li??(i,j?0,1,...,n)? (4.11)

?1,j?i?再令

pn(x)?Ln(x)?y0l0(x)?y1l1(x)???ynln(x)??yi?i?0j?0j?innx?xjxi?xj (4.12)

Ln(x)称为n次拉格朗日插值多项式。

例4.4 已x0=3, x1=-1, x2=4, f(x0)=5, f(x1)=-3, f(x2)=2，求L2(x)。 解：

由式4.10得

l0(x)?l1(x)?l2(x)?(x?x1)(x?x2)(x?1)(x?4)(x?1)(x?4)???(x0?x1)(x0?x2)(3?1)(3?4)4(x?x0)(x?x2)(x?3)(x?4)(x?3)(x?4)??

(x1?x0)(x1?x2)(?1?3)(?1?4)20(x?x0)(x?x1)(x?3)(x?1)(x?3)(x?1)??(x2?x0)(x2?x1)(4?3)(4?1)5由2次拉格朗日插值多项式4.12得

L2(x)?l0(x)y0?l1(x)y1?l2(x)y2(x?1)(x?4)(x?3)(x?4)(x?3)(x?1)?5??(?3)??2

4205??x2?4x?2??例4.5 根据拉格朗日插值公式，ln(x)在x=1.5的值，差值区间x∈[1,20]。 Matlab int-lagrange

4.4 插值与拟合的其他方法

1 差商与牛顿插值

设在x0,x1,…,xn处，函数y=f(x)的取值分别为f(x0), f(x1),…, f(xn)。根据pn(x)=a0+a1(x- x0)+…+an(x- x0)(x- x1)…(x- xn-1)取： Nn(x0)= a0=f(x0)

Nn(x1)= a0+a1(x1- x0)=f(x1)

Nn(x2)= a0+a1(x2- x0)+ a2(x2- x0) (x2- x1)=f(x2)

Nn(xn)= a0+a1(xn- x0)+…+ an(xn- x0)(xn- x1)…(xn- xn-1)=f(xn) (4.13) 这是关于未知数a0,a1,…,an的下三角形方程组，可以求得

a0?f(x0)a1?a2?f(x1)?f(x0)?f[x0,x1]x1?x0f(x2)?f(x0)?a1(x2?x0)f[x0,x2]?f[x0,x1]??f[x0,x1,x2](x2?x0)(x2?x1)x2?x1

利用数学归纳法可以证明ak?f[x0,x1,?,xk](k?1,2,...,n)，代入式4.13便得到

Nn(x)?f(x0)?f[x0,x1](x?x0)?f[x0,x1,x2](x?x0)(x?x1)?...?f[x0,x1,?,xn](x?x0)(x?x1)...(x?xn?1)式4.14称为n次牛顿(Newton)插值公式。 例4.6 利用牛顿插值法求解例4.4。 解：

(4.14)

f[x0]?f(x0)?5f[x0,x1]?f[x1,x2]?f(x1)?f(x0)?3?5??2

x1?x0?1?3f(x2)?f(x1)2?3??1x2?x14?1f[x1,x2]?f[x0,x1]1?2???1x2?x04?3根据差商的性质，有

f[x0,x1,x2]?N2(x)?f(x0)?f[x0,x1](x?x0)?f[x0,x1,x2](x?x0)(x?x1)

?5?2(x?3)?(?1)(x?3)(x?1)??x2?4x?22 列维尔插值法

列维尔(Neville)插值法又称为逐次插值，是一种实用的插值方法。其基本思想是通过低一次的多项式的组合来得到高一次的插值多项式。 列维尔插值法的算法顺序如下。

P0,1是通过比它低一次的插值多项式P0和P1与(x-x0)和(x-x1)交叉得到的，即

P0,1?(x?x0)Px?x0x?x11?(x?x1)P0?P?P0 1(x?x0)?(x?x1)x1?x0x0?x1同样，

P1,2?(x?x1)P2?(x?x2)Px?x2x?x11?P?P2 1(x?x1)?(x?x2)x1?x2x2?x1P0,1,2是P0,1和P1,2与(x-x0)和(x-x2)交叉得到的，即

P0,1,2(x?x0)Px?x2x?x01,2?(x?x2)P0,1??P0,1?P1,2(x?x0)?(x?x2)x0?x2x2?x0(x?x1)P2,3?(x?x3)Px?x3x?x11,2?P?P2,3 1,2(x?x1)?(x?x3)x1?x3x3?x1(x?x0)Px?x3x?x01,3?(x?x3)P0,2?P0,2?P1,3(x?x0)?(x?x3)x0?x3x3?x0P1,2,3?P0,1,2,3?列维尔插值法的插值过程见表4.1。

表4.1列维尔插值法的插值过程 xi x0 x1 f(xi) p0 p1 一阶多项式 p0,1 二阶多项式 三阶多项式 n阶多项式 x2 x3 . . . xn p2 p3 . . . pn p1,2 p2,3 . . . pn-1,n p0,1,2 p1,2,3 . . . pn-2,n-1,n p0,1,2,3 . . . pn-3,n-2,n-1,n p0,1,2,…,n 例4.7 考虑第一类零阶贝塞尔(Bessel)函数J0，下表列出了J0在几个点上的函数值： xi f(xi) 1.0 0.7651977 1.3 0.6200860 1.6 0.4554022 1.9 0.2818186 2.2 0.1103623 2.5 -0.0483838 试用列维尔插值法和拉格朗日插值法计算f(1.5)的近似值。

3 曲线拟合的最小二乘法

最小二乘法应用广泛，其原理是建立某种误差的平方和，通过使误差平方和最小，求出问题的解。

例4.8 已知一超定方程组：

2x1?4x2?13x1?5x2?3x1?2x2?62x1?x2?7用最小二乘法进行求解。 解：

建立下面的误差方程：

Q?(2x1?4x2?1)2?(3x1?5x2?3)2?(x1?2x2?6)2?(2x1?x2?7)2

通过将误差平方和Q对x1，x2求一阶导数，并令其为0，得到两个关于x1，x2的线性方程，通过求解该线性方程组，可得到超定方程的解。 最小二乘法数据拟合为：对于一组观测数据（xi,yi）（i=0，1，… ,n），要求出自变量x与因变量y的函数关系y=p(x,a0, a1,…, am)，其中ak（k=0,1,2,… ,m）为待定参数，使给定点xi上的误差?i?f(xi)?yi的平方和

??i?0n2i最小。进求一多项式

p(x)?a0?0(x)?a1?1(x)???am?m(x),m?n (4.15)

使

I(a0,a1,...,am)???i[p(xi)?yi]2i?0n???i[a0?0(xi)?a1?1(xi)?...?am?m(xi)?yi]2i?0n

为最小。这就是最小二乘法逼近，得到的拟合曲线，这种方法称为曲线拟合的最小二乘法。由极值必要条件可得

n?I?2??i[a0?0(xi)?a1?1(xi)?...?am?m(xi)?yi]?k(xi)?0 (4.17) ?aki?0k?0,1,...,m

本文来源：https://www.bwwdw.com/article/3mn6.html