清华大学数学实验_实验9 非线性规划1
高教数学建模与数学实验第 非线性规划PPT课件
(3)
以
X
k 1
D0为初始点，求解min X D 0
I
X
,
rk
，其中
X
D0的
最优解设为 X k = X rk D0；
(4)
检验是否满足
r
m
ln
i=1
gi
Xk
或
rk
m
i=1gi
1
X
，若满
足，停止迭代，令 X * X k ；否则取rk1 = rk ，令k = k 1，
返回
罚函数法
罚函数法基本思想是通过构造罚函数把 约束问题转化为一系列无约束最优化问题， 进而用无约束最优化方法去求解．这类方法 称为序列无约束最小化方法．简称为SUMT 法．
其一为SUMT外点法，其二为SUMT内点
法．
第5页/共43页
SUTM外点法
对一般的非线性规划： min f X
s.t.
k j
j = 1,L, n
求解该线性规划问题，得到最优解X k1 ；
(4) 检验 X k1对原约束是否可行.若 X k1对原约束可行，则转
步骤(5)；否则，缩小步长限制，令
k j
=
k j
j =1,L,n，返
回步骤(3)，重解当前的线性规划问题；
5)
判断精度：若
k j
j =1,L,n，则点 X k1为近似最优解；
这里的罚函数只对不满足约束条件的点实行惩罚:当X D 时，满足
各 gi X 0, hi X = 0 ，故罚项为0，不受惩罚．当X D 时,必
有约束条件 gi X 0或hi X 0 ，故罚项大于0，要受惩罚．
非线性规划
非线性规划非线性规划是一种涉及非线性目标函数和/或非线性约束条件的优化问题。
与线性规划不同，非线性规划可能存在多个局部最优解，而不是全局最优解。
非线性规划在许多领域都有广泛的应用，如经济学、工程学和管理学等。
非线性规划的一般形式可以表示为：最小化或最大化 f(x)，其中 f(x) 是一个非线性函数，x 是决策变量向量。
满足一组约束条件g(x) ≤ 0 和 h(x) = 0，其中 g(x) 和 h(x) 是非线性函数。
为了求解非线性规划问题，可以使用不同的优化算法，如梯度下降法、牛顿法和拟牛顿法等。
这些算法的目标是找到目标函数的最小值或最大值，并满足约束条件。
非线性规划的难点在于寻找全局最优解。
由于非线性函数的复杂性，这些问题通常很难解析地求解。
因此，常常使用迭代算法来逼近最优解。
非线性规划的一个重要应用是在经济学中的生产计划问题。
生产活动通常受到多个因素的限制，如生产能力、原材料和劳动力等。
非线性规划可以帮助确定最佳的生产数量，以最大化利润或最小化成本。
另一个应用是在工程学中的优化设计问题。
例如，优化某个结构的形状、尺寸和材料以满足一组要求。
非线性规划可以帮助找到最佳设计方案，以最大程度地提高性能。
在管理学中，非线性规划可以用于资源分配和风险管理问题。
例如，优化一个公司的广告预算，以最大程度地提高销售额。
非线性规划可以考虑多种因素，如广告投入和市场需求，以找到最佳的广告投放策略。
总之，非线性规划是一种重要的优化方法，用于解决涉及非线性目标函数和约束条件的问题。
它在经济学、工程学和管理学等领域有广泛的应用。
尽管非线性规划的求解难度较大，但通过合适的优化算法，可以找到最佳的解决方案。
第十讲非线性规划一运筹学清华大学林谦
·凸函数：定义在凸集上的函数f(X)称为凸函数，条件是 对于每一对x1，x2及每一个a，0≤a≤1存在：
f(ax1+(1－a)x2)≤a f(x1)+1(1－a)f(x2)
上式中，若≤变为<，则称为严格凸函数。
page 12 3 August 2019
Prof. Wang School of Economics & Management
B) 对于所有d，则dT▽2 f(X*)·d≥0
ii)判断极小点的充分条件
命题（二阶充分条件——无约束）：设f(X)C2 是定义在 以X*为内点的一个区域上的函数，若
A) ▽f(X*)=0 B) Hess阵H(X*)正定（或负定）
则X*是f(X)的严格极小点（或极大点）
page 11 3 August 2019
目标函数 约束条件
page 3 3 August 2019
max：f(X) =30x1+450x2
0.5x1+2x2+0.25x22≤800
x1≥0，x2≥0
Prof. Wang School of Economics & Management
Operations Research
第二十一讲
§2 非线性规划的数学模型及 特点 （1）
page 15 3 August 2019
Prof. Wang School of Economics & Management
Operations Research
第二十一讲
§4 非线性规划求解方法分类（1）
1．一维最优化 ①斐波那契（Fibonacci）法 ②黄金分割法（0.618法） ③牛顿法（切线法） ④抛物线逼近法 ⑤成功和失败法
实验9 非线性规划
数学实验（非线性规划）郭明钊 2012011880 化21一、 生产安排问题1、 问题分析：根据题目的要求可以知道这个问题的目标是使公司获得的利润最大，要做的决策就是如何安排生产，约束的条有：原料供应量，市场的需求量，产品的含硫量。
根据这些条件列出约束方程，运用LINGO 就可以很方便的解出最优解。
2、 建立模型：用(11,2)ij x i j ==，2，3；表示第i 种原料（这里1,2,3分别代表甲，乙，丙）用于生产第j 种产品（这里1,2分别代表A ，B 两种产品）。
目标函数：设公司获利为Z 千元，则目标函数为：33232121231111191561610i i j j j i i j j j MaxZ x x x x x ======+---∑∑∑∑∑（总利润=总的售价-总的成本）约束条件：甲、乙、丙三种原料的供应量都不超过500t ，所以有)3,2,1(50021=≤∑=i xj ijA 、B 两种产品的产量分别不能超过100t ，200t 。
所以有311100i i x=≤∑和321200i i x =≤∑ 产品A ，B 的含硫量分别不能超过2.5%，1.5%，且已知甲、乙、丙三种原料的含硫量分别是3%，1%，2%，于是就有1121311121313%1%2% 2.5%x x x x x x ++≤++ （1） 1222321222323%1%2% 1.5%x x x x x x ++≤++ （2） 由（1）（2）两式进一步化简可得1121310.5 1.50.50,x x x --≤和1222321.50.50.50x x x -+≤另外还有，按照生产工艺要求，原料甲、乙必须首先倒入混合池中混合，混合后的液体再分别与原料丙混合生产A ，B 。
所以在A ，B 两种产品中甲，乙两种原料的相对含量是相同的。
于是就有11211222::x x x x =3、 LINGO 求解（1） 对于第一问，直接根据上面模型在lingo 中输入以下内容：max=9*(x11+x21+x31)+15*(x12+x22+x32)-6*(x11+x12)-16*(x21+x22)-10*(x31+x32); x11+x12<=500;x21+x22<=500;x31+x32<=500;x11+x21+x31<=100;x12+x22+x32<=200;x11*x22-x21*x12=0;0.5*x11-1.5*x21-0.5*x31<=0;1.5*x12-0.5*x22+0.5*x32<=0;得到如下结果：Global optimal solution found.Objective value: 400.0000Objective bound: 400.0000Infeasibilities: 0.000000Extended solver steps: 2Total solver iterations: 131Model Class: NLPTotal variables: 6Nonlinear variables: 4Integer variables: 0Total constraints: 9Nonlinear constraints: 1Total nonzeros: 28Nonlinear nonzeros: 4Variable Value Reduced Cost X11 0.000000 0.000000X21 0.000000 4.000000X31 0.000000 0.000000X12 0.000000 2.000000X22 100.0000 0.000000X32 100.0000 0.000000Row Slack or Surplus Dual Price1 400.0000 1.0000002 500.0000 0.0000003 400.0000 0.0000004 400.0000 0.0000005 100.0000 0.0000006 0.000000 2.0000007 0.000000 0.2000000E-018 0.000000 2.0000009 0.000000 6.000000于是可以得出结论，在第一问的条件下，最优的生产安排应该是用100t的乙原料和100t的丙原料混合，生成200t产品B，所得的最大利润为40万元（400千元）。
实验九(非线性规划)
大学数学实验实验报告——非线性规划2 大学数学实验 实验报告 | 2014/5/22x3A +x3B ≤500;xA ≤100; xB ≤200;x1A ,x2A ,x3A ,x1B ,x2B ,x3B ,xA ,xB ≥0;针对于不同的情况，改变约束条件中的相关参数，就可以进行优化求解了。
解决方案：针对第一问，直接按照上面的不等式利用LINGO 软件直接编写程序求解： max =-6*x1A-16*x2A-10*x3A-6*x1B-16*x2B-10*x3B+9*xA+15*xB; x1A+x2A+x3A=xA; x1B+x2B+x3B=xB;3*x1A+x2A+2*x3A<=2.5*xA; 3*x1B+x2B+2*x3B<=1.5*xB; x1A*x2B=x2A*x1B; x1A+x1B<=500; x2A+x2B<=500; x3A+x3B<=500; xA<=100; xB<=200; 非负约束省略。
运行得到优化结果报告的部分内容如下： Global optimal solution found.Objective value: 400.0000 Objective bound: 400.0000 Infeasibilities: 0.000000 Extended solver steps: 2 Total solver iterations: 134Model Class: NLPVariable Value Reduced Cost X1A 0.000000 0.000000 X2A 0.000000 4.000000 X3A 0.000000 0.000000 X1B 0.000000 2.0000003大学数学实验 实验报告 | 2014/5/2X2B 100.0000 0.000000 X3B 100.0000 0.000000 XA 0.000000 0.000000 XB 200.0000 0.000000Row Slack or Surplus Dual Price 1 400.0000 1.000000 2 0.000000 -14.00000 3 0.000000 -22.00000 4 0.000000 2.000000 5 0.000000 6.000000 6 0.000000 0.2000000E-01 7 500.0000 0.000000 8 400.0000 0.000000 9 400.0000 0.000000 10 100.0000 0.000000 11 0.000000 2.000000从这份报告可以看出，当用乙、丙各100t 来生产200t 的B 时利润最大为40（万元）。
数学实验报告-非线性规划与多目标规划实验
1）建立函数M文件：
functionf=fun(x)
f=-20*exp(-0.2*(0.5*(x(1)^2+x(2)^2)^0.5))-exp(0.5*(cos(2*pi*x(1))+cos(2*pi*x(2))))+22.713;
2）
x0=[0,0];
options=optimset('display', 'iter', 'tolfun',1e-10);
1982
0.117
0.465
0.215
0.187
0.213
0.311
-0.019
0.084
1983
0.092
-0.015
0.224
0.235
0.217
0.08
0.237
-0.128
1984
0.103
0.159
0.061
0.03
-0.097
0.15
0.074
-0.175
1985
0.08
0.366
0.316
[3]熟悉MATLAB软件求解非线性规划模型的基本命令；
[4]通过范例学习，了解建立非线性规划模型的全过程，与线性规划比较其难点何在。
本实验包括基础实验、应用实验和创新实验，基础实验和应用实验要求独立完成，创新实验要求合作完成。通过该实验的学习，使学生掌握最优化技术，认识面对什么样的实际问题，提出假设和建立优化模型，并且使学生学会使用MATLAB软件和Lingo软件求解非线性规划模型，注意初始解的选择不同会导致软件求出的解的变化（是局部最优解还是整体最优解）。解决现实生活中的最优化问题是本科生学习阶段中一门重要的课程，因此，本实验对学生的学习尤为重要。
非线性规划
非线性规划什么是非线性规划？非线性规划（Nonlinear Programming，简称NLP）是一种数学优化方法，用于求解包含非线性约束条件的优化问题。
与线性规划不同，非线性规划中的目标函数和约束条件都可以是非线性的。
非线性规划的数学表达式一般来说，非线性规划可以表示为以下数学模型：minimize f(x)subject to g_i(x) <= 0, i = 1, 2, ..., mh_j(x) = 0, j = 1, 2, ..., px ∈ R^n其中，f(x)是目标函数，g_i(x)和h_j(x)分别是m个不等式约束和p个等式约束，x是优化变量，属于n维实数空间。
非线性规划的解法由于非线性规划问题比线性规划问题更为复杂，因此解决非线性规划问题的方法也更多样。
以下列举了几种常用的非线性规划求解方法：1. 数值方法数值方法是最常用的非线性规划求解方法之一。
它基于迭代的思想，通过不断优化目标函数的近似解来逼近问题的最优解。
常见的数值方法有梯度下降法、牛顿法、拟牛顿法等。
2. 优化软件优化软件是一类针对非线性规划问题开发的专用软件，它集成了各种求解算法和优化工具，可以方便地求解各种类型的非线性规划问题。
常见的优化软件有MATLAB、GAMS、AMPL等。
3. 线性化方法线性化方法是一种将非线性规划问题转化为等价的线性规划问题的求解方法。
它通过线性化目标函数和约束条件，将非线性规划问题转化为线性规划问题，然后利用线性规划的求解方法求解得到最优解。
4. 分类方法分类方法是一种将非线性规划问题分解为若干个子问题求解的方法。
它将原始的非线性规划问题分解为多个子问题，然后将每个子问题分别求解，并逐步逼近原始问题的最优解。
以上仅是非线性规划求解方法的一小部分，实际上还有很多其他的方法和技巧可供选择。
在实际应用中，选择合适的方法和工具是非常重要的。
非线性规划的应用非线性规划在实际生活和工程中有着广泛的应用。
《非线性规划》课件
非线性规划的约束条件
非线性规划的约束条件是指限制问题解的一组方程或不等式。这些约束条件可以包括物理限制、资源约 束和行为限制等。
非线性规划的求解方法
线性化方法
将非线性问题转化为等价的 线性问题，然后使用线性规 划方法求解。
牛顿法
使用牛顿迭代法逐步逼近最 优解。
拟牛顿法
使用近似Hessian矩阵的方法 优化牛顿法。
变尺度法、全局优化方法
1
变尺度法
通过改变尺度，将问题转化为更易求解的形式。
2
全局优化方法
使用启发式算法寻找全局最优解。
非线性规划的应用领域
生产计划问题
优化生产计划，提高效率和利润。
交通运输问题
优化交通网络和运输流程。
优化电力系统
使电力系统运行更加高效和可靠。
决策支持系统
为决策者提供优化建议和决策支持。
医资源分配和治疗方案。
非线性规划的挑战
复杂的问题结构和求解困难。
未来的研究方向
未来的研究方向包括改进算法性能、适用于大规模问题的方法和考虑不确定性的优化模型等。
《非线性规划》PPT课件
在这个《非线性规划》PPT课件中，我们将深入探讨非线性规划的各个方面， 并介绍其在不同领域的应用。让我们一起开启这个激动人心的学习之旅！
什么是非线性规划？
非线性规划是一种在优化问题中寻找最优解的数学方法。它处理的是有非线 性约束条件和目标函数的优化问题。
非线性规划的优化目标
数学实验非线性规划.docx
《大学数学实验》作业非线性规划班级:姓名:学号: 日期:目录【实验目的】 (3)【实验内容】 (3)题目1 （课本习题第九章第4题） (3)【第（1）问求解】 (3)【第（2）问求解】 (7)【第（3）问求解】 (7)【拓展实验、思考、对比、分析】 (8)【木题小结】 (10)题目2（课本习题第九章第8题） (10)【模型建立】 (11)【模型求解】 (14)【第（1）问求解】 (14)【第（2）问求解】 (20)【第（3）问求解】 (22)【拓展实验、思考、对比、分析】 (23)【本题小结】 (25)【实验心得、体会】 (25)注：本实验作业脚本文件均以ex9_4_l形式命名，其中ex代表作业，9_4_1表示第九章第四小题第一个程序。
自编函数均以exf9_4_l形式命名，exf代表作业函数，9_4_1 表示第九章第四题第一个自编函数。
【实验目的】1.掌握用MATLAB优化工具箱和LINGO解非线性规划的方法；2.练习建立实际问题的非线性规划模型。
【实验内容】题目1 （课本习题第九章第4题）某公司将3种不同含硫量的液体原料（分别记为甲、乙、丙）混合生产两种产品（分别记为A, B）。
按照生产工艺的要求，原料甲、乙必须首先倒入混合池中混合，混合后的液体再分别与原料丙混合牛产A, Bo已知原料甲、乙、丙的含硫量分别是3%, 1%, 2%,进货价格分别为6千元/t, 16千元/t, 10千元/t；产品A, B的含硫量分别不能超过2.5%, 1.5%,售价分别为9千元/t, 15千元/t。
根据市场信息，原料甲、乙、丙的供应量都不能超过500t；产品A, B的最大市场需求量分别为100t, 200to（1）应如何安排生产？⑵如果产品A的最大市场需求量增长为600t,应如何安排生产？⑶如果乙的进货价格下降为13千元/t,应如何安排生产？分别对（1）、（2）两种情况进行讨论。
【第（1）问求解】【模型建立】⑴模型该题为带约束非线性规划问题，其模型包含决策变量、FI标函数和约束条件。
非线性实验报告
非线性实验报告非线性实验报告摘要：本实验旨在研究非线性系统的特性，并通过实验验证非线性系统的存在和影响。
实验过程中，我们采用了不同的实验方法和工具，包括数学模型、实验仪器和数据分析软件。
通过实验结果的分析和对比，我们得出了一些关于非线性系统的结论，并对实验中可能存在的误差和限制进行了讨论。
引言：非线性系统是指其输入与输出之间的关系不符合线性关系的系统。
在现实世界中，非线性系统无处不在，如生物系统、电子电路、经济系统等。
了解和研究非线性系统的特性对于我们理解和应用这些系统具有重要意义。
本实验旨在通过实际操作和数据分析，探索非线性系统的行为和特性。
实验方法：我们选择了一种简单的非线性系统作为研究对象，即二次函数。
通过调整二次函数的系数和参数，我们可以观察到不同的非线性行为。
在实验中，我们使用了一台计算机和数据采集卡作为实验仪器，利用数学建模和数据分析软件进行数据处理。
实验步骤：1. 设计二次函数模型：我们首先根据实验要求设计了一个二次函数模型，包括系数和参数的选择。
这个模型可以模拟实际系统中的非线性行为。
2. 数据采集：我们通过计算机和数据采集卡采集了一系列输入和输出数据。
输入数据是实验中施加在系统上的不同信号，输出数据是系统对这些信号的响应。
3. 数据处理和分析：我们使用数据分析软件对采集到的数据进行处理和分析。
首先，我们绘制了输入-输出曲线，以观察系统的非线性特性。
然后，我们对数据进行了拟合和回归分析，以确定二次函数的系数和参数。
实验结果：通过实验和数据分析，我们得到了以下结果：1. 非线性特性的存在：我们观察到系统的输入-输出曲线不是一条直线，而是呈现出弯曲的形状。
这表明系统存在非线性特性。
2. 参数对系统行为的影响：我们发现，调整二次函数的系数和参数可以改变系统的响应。
例如，增加二次项的系数可以使曲线更加陡峭，而增加线性项的系数可以使曲线更加平缓。
3. 非线性现象的局限性：我们也观察到，在一定范围内，系统的响应是线性的。
第十讲--非线性规划(一)(运筹学基础-清华大学-王永县)解析PPT课件
page 9 4/24/2020
Prof. Wang School of Economics & Management
Operations Research
第二十）：设是En子集，f(X) C1(C1表 明存在一阶导数)是上函数，若X*是f( X)在上一个相对
Operations Research
第二十一讲
§2 非线性规划的数学模型及特 点 （3）
[例4-4]非线性规划为
min f(X)=(x1－2)2+(x2－2)2 h(X)=x1+x2－6≤0
显然，此时的最优解为C点(x1*=x2*=2 ，f(X*)=0)，该点落在可 行或内部，其边界约束失去作用。
极小点，则对于任一X*的可行方向dEn必有▽f(X*)·d≥0。 （其中，▽f(X*)表示函数f( X)的一阶梯度或导数）
命题2 （二阶必要条件——有约束情况）：设是En的一 个子集，且f( X) Ｃ2（Ｃ2表明存在二阶导数）是上的一 个函数。若X*是f( X)在上的一个相对极小点。则对于任 一X*处的可行方向dＥn有：
x1≥0，x2≥0
Prof. Wang School of Economics & Management
Operations Research
第二十一讲
§2 非线性规划的数学模型及特 点 （1）
非线性规划的数学模型通常表示成以下形式。
min f(X)
hi(X)=0 i=1，2，…，m gj(X)≥0 j=1，2，…，l
求：公司为获取最大营业额（销售额）的最优营业计划。
page 2 4/24/2020
Prof. Wang School of Economics & Management
数学建模—非线性规划实验报告
实验六数学建模—非线性规划实验目的：1．直观了解非线性规划的基本内容．2．掌握用数学软件求解优化问题．实验内容：1、某厂向用户提供发动机，合同规定，第一、二、三季度末分别交货40台、60台、80台．每季度的生产费用为()2bxaxxf+=（单位:元）, 其中x是该季度生产的台数．若交货后有剩余，可用于下季度交货，但需支付存储费，每台每季度c元．已知工厂每季度最大生产能力为100台，第一季度开始时无存货，设a=50、b=0.2、c=4，问:工厂应如何安排生产计划，才能既满足合同又使总费用最低．讨论a、b、c变化对计划的影响，并作出合理的解释．2、一基金管理人的工作是: 每天将现有的美元、英镑、马克和日元四种货币按当天汇率相互兑换，使在满足需要的条件下，按美元计算的价值最高．设某天的汇率、现有货币和当天需求如下：问该天基金管理人应如何操作. (“按美元计算的价值”指兑入、兑出汇率的平均值，如1英镑相当于()258928.01697.1+=1.696993美元.)实验过程与结果：1、（1）模型建立决策变量：设第1，2，3季度分别生产x1，x2，x3台发动机，第1，2季度末分别有存货40-x1，x1+x2-100台，第3季度末无存货目标函数：设总费用为z=a(x1+x2+x3)+b(x1^2+x2^2+x3^2)+c[(x1-40)+(x1+x2-100)]约束条件：生产的发动机应该在第3季度末全部卖出，则有x1+x2+x3=180；同时要保证第1，2季度能供货且有能力生产，要求x1≥40，x1+x2≥100,100≥x1,100≥x2,100≥x3非负约束：x1,x2,x3≥0综上可得：Maxz=a(x1+x2+x3)+b(x1^2+x2^2+x3^2)+c[(x1-40)+(x1+x2-100)]s.t．x1+x2+x3=180x1+x2≥100x1≥400≤x1,x2,x3≤100（2）模型求解结果为：即工厂应第一季度生产50台发动机，第二季度生产60台发动机，第三季度生产70台发动机，才能既满足合同又使总费用最低。
数学实验报告——利用MALTAB进行非线性规划
㈡简要分析
本题是一道比较简单的非线性规划求解问题， 不涉及数学建模以及应用分析 等问题，只需要将所给问题转化为 MATLAB 代码进行处理。这里可以采用不同的 算法，并对比进行分析。
㈢方法与公式
1、求解规划方法
求解本题可以考虑几种规划。 对于第(1)组约束来说，仅含有上下界约束，可以考虑使用 SQP 方法或者置 信域方法；但是考虑到本题规模较小，这里仅采用 SQP 方法。（事实是，设定使 用大规模算法后，MATLAB 仍旧自动使用 SQP 完成计算，因而无法进行对比）； 对第(2)组约束来说，仅可以使用 SQP 方法； 对第(3)组约束来说，可以使用 SQP 方法，也可以利用后两个等式把规划化 化简后再使用 SQP 方法。
5、脚本 2
5
v1 = -10*ones(1,4); v2 = 10*ones(1,4); x00 = [-3,-1,-3,-1]; x01 = -[-3,-1,-3,-1]; x02 = [10,-20,5,3]; x03 = [-5,10,15,-20]; opt = optimset('largeScale','off','MaxFunEvals',4000,'MaxIter', 1000,'Algorithm','active-set'); [x(1,:),f1,exitflag1,out(1)] = fmincon('fun',x00,[],[],[],[],v1,v2,@edge3,opt); [x(2,1:2),f2,exitflag2,out(2)] = fmincon('fun1',x00(2:3),[],[],[],[],v1(1:2),v2(1:2),@edge 31,opt); [x(3,:),f3,exitflag3,out(3)] = fmincon('fun',x01,[],[],[],[],v1,v2,@edge3,opt); [x(4,1:2),f4,exitflag4,out(4)] = fmincon('fun1',x01(2:3),[],[],[],[],v1(1:2),v2(1:2),@edge 31,opt); [x(5,:),f5,exitflag5,out(5)] = fmincon('fun',x02,[],[],[],[],v1,v2,@edge3,opt); [x(6,1:2),f6,exitflag6,out(6)] = fmincon('fun1',x02(2:3),[],[],[],[],v1(1:2),v2(1:2),@edge 31,opt); [x(7,:),f7,exitflag7,out(7)] = fmincon('fun',x03,[],[],[],[],v1,v2,@edge3,opt); [x(8,1:2),f8,exitflag8,out(8)] = fmincon('fun1',x03(2:3),[],[],[],[],v1(1:2),v2(1:2),@edge 31,opt); answer = zeros(1:8); for i = 1:8 out(i) end for j = 1:4 x(2*j,3)=x(2*j,1); x(2*j,1) = -x(2*j,2);
实验9非线性规划
实验9非线性规划实验目的1.掌握用MATLAB优化工具箱和LINGO解非线性规划的方法；2.练习建立实际问题的非线性规划模型。
实验内容1.对问题增加以下条件，并分别取初值(-3,-1,-3,-1)和(3,1,3,1)，求解非线性规划：(1);(2);(3)再试取不同的初值或用分析梯度计算，比较计算结果。
你能从中得到什么启示？解对于这个非线性规划问题，先对目标函数编写如下”f.m”文件：function[f,g,H]=f(x)f=100*(x(2)-x(1).^2).^2+(1-x(1)).^2+90*(x(4)-x(3).^2).^2+(1-x(3)) .^2+10.1*((1-x(2)).^2+(1-x(4)).^2)+19.8*(x(2)-1)*(x(4)-1);if nargout>1g(1)=-400*(x(2)-x(1).^2)*x(1)-2*(1-x(1));g(2)=200*(x(2)-x(1).^2)-20.2*(1-x(2))+19.8*(x(4)-1);g(3)=-360*(x(4)-x(3).^2)*x(3)-2*(1-x(3));g(4)=180*(x(4)-x(3).^2)-20.2*(1-x(4))+19.8*(x(2)-1);endif nargout>2H=[1200*x(1).^2-400*x(2)+2 -400*x(1) 0 0-400*x(1) 220.2 0 19.80 0 1080*x(3).^2-360*x(4)+2 -360*x(3)0 19.8 -360*x(3) 200.2];end对于后两个约束条件编写函数”c_2.m”、”c_3.m”：%c_2.mfunction[c1,GC1]=c_2(x)c1=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10;x(1)*x(2)*x(3)*x(4)-100 ;-x(1)*x(2)*x(3)*x(4)-100];if nargout>2GC1=[x(2)-1 x(1)-1 0 0-x(2) -x(1) 0 0x(2)*x(3)*x(4) x(1)*x(3)*x(4) x(2)*x(1)*x(4) x(2)*x(3)*x(1) -x(2)*x(3)*x(4) -x(1)*x(3)*x(4) -x(2)*x(1)*x(4)-x(2)*x(3)*x(1)]’;end%c_3.mfunction[c1,c2,GC1,GC2]=c_3(x)c1=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10];c2=[x(1)+x(2);x(1)*x(2)*x(3)*x(4)-16];if nargout>2GC1=[x(2)-1 x(1)-1 0 0-x(2) -x(1) 0 0]’;GC2=[1 1 0 0x(2)*x(3)*x(4) x(1)*x(3)*x(4) x(2)*x(1)*x(4)x(2)*x(3)*x(1)]’;end(1)取初值(-3,-1,-3,-1)，首先用中等规模算法，分别采用数值梯度和分析梯度计算，程序代码如下：x0=[-3 -1 -3 -1]; %赋初值opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000);%采用数值梯度计算[x,fv,ef,out,lag,grad,hess]=fmincon(@f,x0,[],[],[],[],-10*ones(1,4),10*ones(1,4),[],opt1)pauseopt2=optimset(opt1,'GradObj','on','GradCon','on','DerivativeCheck ','on'); %采用分析梯度计算[x,fv,ef,out,lag,grad,hess]=fmincon(@f,x0,[],[],[],[],-10*ones(1,4),10*ones(1,4),[],opt2)由于约束中仅含有上下界约束，所以可以使用大规模算法：x0=[-3 -1 -3 -1]; %赋初值opt1=optimset('MaxIter',3000,'MaxFun',20000); %采用数值梯度计算[x,fv,ef,out,lag,grad,hess]=fmincon(@f,x0,[],[],[],[],-10*ones(1,4),10*ones(1,4),[],opt1)pauseopt2=optimset(opt1,'GradObj','on','GradCon','on','DerivativeCheck ','on'); %采用分析梯度计算[x,fv,ef,out,lag,grad,hess]=fmincon(@f,x0,[],[],[],[],-10*ones(1,4),10*ones(1,4),[],opt2)在MATLAB上分别运行上述两个程序，即可以得到分别采用大规模算法、中等规模算法以及分析梯度、数值梯度所得到的结果，将各种算法所得的值、得到的目标函数最优值、目标函数调用次数放在一起中作对比，可以得到如下表最优值（需要说明的是，前三种算法得到的的值并不完全相同，只是在保留6位小数之后看起来一样。
非线性规划
非线性规划非线性规划(Nonlinear Programming ，简记为NP)研究的对象是非线性函数的数值最优化问题，是运筹学的最重要分支之一，20世纪50年代形成一门学科，其理论和应用发展十分迅猛，随着计算机的发展，非线性规划应用越来越广泛，针对不同的问题提出了特别的算法，到目前为止还没有适合于各种非线性规划问题的一般算法，有待人们进一步研究.§1 非线性规划基本概念一、引例例7.1 一容器由圆锥面和圆柱面围成. 表面积为S ，圆锥部分高为h ，h 和圆柱部分高2x 之比为a ,1x 为圆柱底圆半径.求21,x x 使面积最大.解： 由条件 a x h =2/22121231x x x ax V ππ+=21212222112221x x x x a x x S πππ+++⋅⋅=所以，数学模型为：212)311(max x x a V π+=s.t. S x x x x a x x =+++21212222112πππ0,21≥x x例7.2 某高校学生食堂用餐，拟购三种食品，馒头0.3元/个，肉丸子1元/个，青菜0.6/碗.该学生的一顿饭支出不能够超过5元.问如何花费达到最满意？解： 设该学生买入馒头，肉丸子，青菜的数量分别为321,,x x x ； 个人的满意度函数即为效用函数为321321321),,(aaax x Ax x x x u =.于是数学模型为321321321),,(max aaax x Ax x x x u =s.t.56.03.0321≤++x x x 321,,x x x 为非负整数二、数学模型称如下形式的数学模型为数学规划（Mathematical Programming 简称MP ） )(min x f z = （7.1） （MP ） t s . 0)(≥x g i m i ,,1 = （7.2） 0)(=x h j l j ,,1 = （7.3）其中Tn x x x x ),,,(21 =是n 维欧几里得空间nR 中的向量（点），)(x f 为目标函数，0)(,0)(=≥x h x g j i 为约束条件.称满足约束条件的向量x 为（MP ）问题的一个可行解，全体可行点组成的集合称为可行域.K ={}l j x h mi x g R x j i n,,2,10)(,,2,10)( ===≤∈如果)(),(),(x h x g x f j i 均为线性函数，就是前面所学的线性规划问题(LP).如果至少有一个为非线性函数称为非线性规划问题.由于等式约束0)(=x h j 等价于下列两个不等式约束 0)(,0)(≥-≥x h x h j j 所以(MP)问题又可表示为 )(min x f z =s.t. 0)(≥x g i m i ,,1 = （7.4） 三、数学基础 1、线性代数知识考虑二次型Az z T ，z 为n 维向量正定的二次型：若对于任意0≠z ，有0>Az z T； 半正定的二次型：若对于任意0≠z ，有0≥Az z T ； 负定的二次型：若对于任意0≠z ，有0<Az z T ； 半负定的二次型：若对于任意0≠z ，有0≤Az z T ；不定二次型：0≠∃z ，有0>Az z T，又0≠∃z ，有0<Az z T.二次型Az z T 为正定的充要条件是它的矩阵A 的左上角各阶主子式都大于零. 二次型Az z T 为负定的充要条件是它的矩阵A 的左上角各阶主子式负正相间.2、分析数学知识（1）方向导数和梯度(二维为例)考虑函数),(21x x f Z =，及方向j i lϕϕsin cos +=梯度：Tx f x f j x f i x f x x f ),(),(212121∂∂∂∂=∂∂+∂∂=∇ ； 方向导数：⎪⎪⎭⎫⎝⎛∂∂∂∂=∂∂+∂∂=∂∂ϕϕϕϕsin cos ),(sin cos 2121x f x f x f x f l f )),,(cos(||),(||),(),(21212121l x x gardf x x gardf lx x gardf lx x f T=⋅=⋅∇=考虑等值线00201),(c x x f =上一点),(0201x x 梯度方向 ),(0201x x gardf 即为法线方向n.如果)(x f 二次可微，称⎪⎪⎪⎪⎪⎭⎫⎝⎛=)()()()()()()()()()(212222111211x h x h x h x h x h x h x h x h x h x H nn n n n n为)(x f 在点 x 处的Hesse 矩阵.（2）多元函数泰勒公式：若)(,),(0x f R S x x f u n⊆∈=在点0x 处的某个领域具有二阶连续偏导数，则有x x x f x x x f x f x x f T T∆∆+∇∆+∆∇+=∆+)(21)()()(02000θ 10≤≤θ )||(||)(21)()(||)(||)()(2020000x x x f x x x f x f x x x f x f T TT ∆+∆∇∆+∆∇+=∆+∆∇+=οο 四、最优解的类型定义7.1 (MP)问题的一个可行点*x 被称为整体极小点，如果对于任意的可行点K x ∈，都有不等式)()(*x f x f ≥成立.如果对于任意可行点*,x x K x ≠∈均有)()(*x f x f >，称点*x 是)(x f 的可行解集K上的严格整体极小点.定义7.2 问题(MP)的一个可行点*x 被称为一个局部极小点，如果存在一个正数ε使得对于所有满足关系式ε<-*x x 的可行点x 都有)()(*≥x f x f 成立.如果对任意的可行点K x ∈，*≠x x ，存在一个正数ε使得对于所有满足关系式ε<-*x x 的可行点x 都有)()(*>x f x f 成立.则称*x 是)(x f 在K 上的一个局部严格极小点.显然整体极小点一定是局部极小点，反之不然. 五、凸规划定义7.3 集合K 被称为nR 中的一个凸集，如果对于K 中任意两点21,x x 和任一实数]1,0[∈λ，点K x x ∈-+21)1(λλ.几何解释：当一个集合是凸集时，连接此集合中任意两点的线段也一定包含在此集合内，规定φ空集是凸集.定义7.4 凸函数：)(x f 是凸集K 上的实值函数，如果对于K 中任意两点21,x x 和任意实数]1,0[∈λ有不等式)()1()())1((2121x f x f x x f λλλλ-+≤-+成立.严格凸函数：)(x f 是凸集K 上的实值函数，如果对于K 中任意两点21,x x ，21x x ≠和任意实数)1,0(∈λ，有不等式)()1()())1((2121x f x f x x f λλλλ-+<-+成立.定义7.5 )(x f 是定义在凸集K 上的实值函数，如果)(x f -是K 上凸函数，称)(x f 是凹函数.定理7.1 设)(x f 是凸集K 上的凸函数，则)(x f 在K 中的任一局部极小点都是它的整体极小点.证明： 设*x 是一局部极小点而非整体极小点，则必存在可行点K x ∈(可行域))()(*x f x f <∍.对任一]1,0[∈λ，由于)(x f 的凸性，有 )()()1()())1((***x f x f x f x x f ≤-+≤-+λλλλ当0→λ时，*)1(x x λλ-+与*x 充分接近，与*x 是局部极小矛盾. 证毕. 定义7.6 设有(MP)问题)(min x f kx ∈，若可行域K 是凸集，)(x f 是K 上的凸函数，则称此规划问题为凸规划.定理7.2 凸规划的任一局部极小解为整体极小解. 六、非线性规划问题的求解方法 考虑(MP)问题：lj x h m i x g t s x f j i ,,10)(,,10)(.)(min ===≥ （7.5） 一般来说，MP 问题难以求得整体极小点，只能求得局部极小点.以后我们说求(MP)问题，指的是求局部极小点.1、无约束极小化问题（UMP ） )(min x f nRx ∈ （7.6） 这里)(x f 是定义在n R 上的一个实值函数定理7.3（一阶必要条件）如果)(x f 是可微函数.*x 是上述无约束问题（UMP ）的一个局部极小点或局部极大点的必要条件是：0)(*=∇x f .满足0)(=∇x f 的点称为平稳点或驻点.定理7.4 设)(x f 为定义在n R 上的二阶连续可微函数，如果*x 是)(x f 的一个局部极小点，必有nT Ry y x H y x f ∈∀≥=∇0)(0)(**这里)(*x H 表示)(x f 在*x 处的Hesse 矩阵.证明：nE y ∈∀，根据)(x f 在点*x 处的展开式有)()(21)()(2*2**λολλ++=+y x H y x f y x f T)0)((*=∇x f若0)(,*<∍∈∃y x H y R y T n ，当λ充分小时，有 )()(21|2*2λολ>y x H y T∴有)()(**x f y x f <+λ.这和*x 是)(x f 的极小矛盾.定理7.5 设)(x f 是定义在nR 上的二阶连续可微函数，如果点*x 满足0)(*=∇x f ，而且存在*x 的一个邻域0)(),(,),(*≥∈∀∈∀∍*y x H y x x R y x T n 有 ，则*x 是)(x f 的一个局部极小点.在高等数学中求解极值点方法先求出满足0)(=∇x f 的点及不可导点.在这些点判断)(x f 是否取得极小值.2、等式约束的极小化问题考虑 )(min x fl j x h t s j ,,10)(. == （7.7）定义7.7 如果)(,),(),(21x h x h x h l ∇∇∇ 在点x 处线性无关，则称点x 是此约束条件的一个正则点.Langrange 乘子法：引进拉格朗日函数 ∑=-=lj jj x h u x f u x L 1)()(),(其中Tl u u u u ),,,(21 =被称为拉格朗日乘子向量.定理7.6 设l j x h x f j ,,1),(),( =是连续可微函数，*x 是)(x f 在可行集中的一个局 部极小点.在*x 是正则点的假定下必存在一个拉格朗日乘子向量u ，使得),(*u x 满足方程组)(0)()(*1**==∇-∇∑=x h x h u x f lj j j对等式约束，用拉格朗日乘子法求解出平稳点，判断是否极值点.用上述解析法求解无约束和等式约束极值问题的平稳点，再判断是否为极值点.该方法有一定的局限性：（1）它们要求函数是连续的，可微的，实际问题中不一定满足这一条件； （2）上述求平稳点的方程组求解比较困难，有些解不出来； （3）实际问题中有大量的不等式约束.因此求解非线性规划应有更好的新方法.实际应用中一般用迭代法来求解非线性规划问题，即求近似最优解的方法.3、非线性规划问题的求解方法—迭代法迭代法一般过程为：在(MP)问题的可行域K 内选择初始点0:,0=k x ，确定某一方向k p ，使目标函数)(x f 从k x 出发，沿k p 方向使目标函数值下降，即)0(,>∈+=λλK p x x k ，有)()(0x f x f <，进一步确定kλ，使)(m i n )(0k k k k k p x f p x f λλλ+=+>，令k k k k p x x λ+=+1.如果1+k x 已满足某终止条件，1+k x 为近似最优解.否则，从1+k x 出发找一个方向1+k p ，确定步长1+k λ，使K p x x k k k k ∈+=++++1112λ，有)(min )(1102++>++=k k k p x f x f λλ.如此继续，将得到点列{}kx .显然有 >>>>)()()(1kx f x f x f ,即点列{}kx 相对应的目标函数是一个单调下降的数列.当{}kx 是有穷点列时，希望最后一个点是(MP)问题的某种意义下的最优解.当{}kx 为无穷点列时，它有极限点，其极限点是(MP)的某种意义下的最优解（此时称该方法是收敛的）.迭代法求解(MP)的注意点：该方法构造的点列{}kx ，其极限点应是近似最优解，即该算法必须是收敛的.迭代法一般步骤：①. 选取初始点0x ，0:=k ②. 构造搜索方向kp ③. 根据kp 方向确定k λ ④. 令k k k k p x xλ+=+1⑤. 若1+k x已满足某终止条件，停止迭代，输出近似最优解1+k x.否则令1:+=k k ，转向第②步.计算终止条件在上述迭代中有：若1+k x满足某终止条件则停止计算，输出近似最优解1+k x.这里满足某终止条件即到达某精确度要求.常用的计算终止条件有以下几个：（1）自变量的改变量充分小时，11||||ε<-+k k x x，或21||||||||ε<-+kk k x x x ，停止计算. （2）当函数值的下降量充分小时，31)()(ε<-+k kx f x f ，或41|)(|)()(ε<-+k k k x f x f x f ， 停止计算.（3）在无约束最优化中，当函数梯度的模充分小时51||)(||ε<∇+k x f ，停止计算.迭代法的关键：① 如何构造每一轮的搜索方向kp ② 确定步长k λ关于k λ，前面是以)(min kk p x f λλ+产生的，也有些算法k λ取为一个固定值，这要根据具体问题来确定.关于kp 的选择，在无约束极值问题中只要是使目标函数值下降的方向就可以了，对于约束极值问题则必需为可行下降方向.定义7.8 设0,,:1≠∈→p R x R R f nn，若存在0>δ使),0(δλ∈∀，)()(x f p x f <+λ则称向量p 是函数)(x f 在点x 处的下降方向.定义7.9 设0,,,≠∈∈∈p R p K x R K nn，若存在0>λ使得K p x ∈+λ，称向量p 是点x 处关于K 的可行方向. 若一个向量p 既是目标函数f 在点x 处的下降方向，又是该点处关于可行域K 的可行方向，则称p 为函数f 在点x 处关于区域K 的可行下降方向.根据不同原理产生了不同的kp 和k λ的选择方法，就产生了各种算法. 在以后的讨论中，一维极值的优化问题是讨论求解步长k λ.无约束极值中讨论的最速下降法，共轭方向法，坐标轮换法，牛顿法，变尺度法及有约束极值中讨论的可行方向法都是根据不同的原理选择kp 得到的迭代算法.七、迭代算法的收敛性定义7.10 设有一算法A ，若对任一初始点（可行点），通过A 进行迭代时，或在有限步后停止得到满足要求的点；或得到一个无穷点列，它的任何一个聚点均是满足要求的点，则称此算法A 具有全局收敛性.定义7.11 设（UMP ）问题的目标函数Px Qx x x f T+=21)(，Q 为对称半正定矩阵， 若由算法A 进行迭代时，不论初始点0x 如何选择，必能在有限步后停止迭代，得到所要求的点，则称此算法A 有二次有限终止性.定义7.12 设序列{}kr收敛于*r，定义满足∞<=--≤**+∞−→−βhkk k rr r r 1______lim0的非负数h 的上确界为{}k r 的收敛级.若序列的收敛级为h ，就称序列是h 级收敛的.若1=h 且1<β就称序列是以收敛比β线性收敛的. 若1>h 或1=h 且0=β就称序列是超线性收敛的. 如何判别算法的收敛性和收敛速度问题本书不讨论.本书给出的算法中，最速下降法具有全局收敛性、是线性收敛的；改进牛顿法具有全局收敛性、二次有限终止性、是二阶线性收敛的；变尺度法和共轭方向法具有全局收敛性和二次有限终止性、是超线性收敛的；Zoutenddijk 法不具有全局收敛性、改进的T-V 法具有全局收敛性；制约函数法具有全局收敛性.关于这些算法的收敛性的讨论和证明有兴趣的读者可参考其他文献.§2 一维极值问题的优化方法前面讨论迭代算法时,从kx 出发确定沿k p 方向的步长k λ是这样求解的),(min 0k k p x f λλ+>这里k x ,k p 已知.所以,若记)()(λλg p x f k k =+,则为求解)(min 0λλg >的过程.于是我们考虑如下形式的极值问题.)(min x f bx a ≤≤ (7.8)b a R x ,,1∈为任意实数很显然可应用高等数学中学过的解析法,即令0)('=x f 求出平稳点,但如前面所述如果该方程解不出来,怎么办?可用下述迭代算法—0.618法.定义7.13 )(x f 定义在],[b a 上,若存在点∍∈],[*b a x 当*x y x ≤<,有)()(y f x f >,当*x y x ≥>时,)()(y f x f >,称)(x f 在],[b a 中为单峰函数(单谷函数).显然满足定义要求的点*x 是)(x f 在],[b a 中的极小点.在],[b a 中任选两点21,x x ,且b x x a <<<21,根据)(x f 的单峰性,若)()(21x f x f <,则*x 必然位于],[2x a 内,如果)()(21x f x f >,则*x 必然位于],[1b x 内.如果)()(21x f x f =,则*x 必然位于],[21x x ,记此区间为],[11b a .如此继续,得闭区间套⊃⊃⊃⊃],[],[],[11n n b a b a b a .显然 ,1,0],,[*=∈i b a x i i ，又0→-i i a b .由闭区间套性质, *x 为极小值点.闭区间套的选择方法不同,求得的*x 的快慢及求解过程的计算量是不一样的，有一个常用的方法是0.618法.0.618法： 取],[],[b a =βα① 在],[βα中选取1λ和2λ,)(618.0),(382.021αβαλαβαλ-+=-+=,求出)(1λf 和)(2λf 进入②.② 若)()(21λλf f <,取],[],[2λαβα=,若αλ-2已足够小,停止,否则进入③.若)()(21λλf f >,取],[],[1βλβα=,若1λβ-已足够小,停止,否则进入④. 若)()(21λλf f =,取],[],[21λλβα=,若12λλ-已足够小,停止,否则进入①. ③ 取上一步中的1λ为2λ,显然有)(618.02αβαλ-+=,令)(382.01αβαλ-+=,求出)(1λf ,返回②.④ 取上一步的2λ为1λ,则有)(382.01αβαλ-+=,令)(618.02αβαλ-+=,求出)(2λf 返回②.设初始区间为],[b a ，用0.618法，经过k 次迭代后],[βα的长度ka b 618.1)(-=-αβ，只要k 充分大αβ-可以小于任何给定的正数.例7.3 用0.618法求解λλλ2)(min 2+=f单峰区间为[-3，5]，2.0=ε解：[α，β]=[-3，5]1λ=-3+0.382×8=0.056 )(1λf =0.1152λ=-3+0.618×8=1.944 )(2λf =7.667由于)(1λf <)(2λf 所以新的不定区间为[α，β] =[-3，1.944] 由于β-α=4.944>0.22λ:=1λ=0.056 )(2λf :=)(1λf =0.115 1λ=-3+0.382×4.944=-1.112 )(1λf =-0.987如此反复得下表7-1：在进行8次迭代后,2.0112.1936.0<+-=-αβ取中间值024.1*-=λ或032.12-=λ作为近似最优解.显然真正极小点是-1.0.一维收索方法很多，如函数逼近法、牛顿法等可参考其他文献.§3 无约束极值的优化方法考虑无约束最优化问题（UMP ）)(min x f nR x ∈ (7.9) 前面已经讨论过,求解此无约束优化问题,可以求出平稳点及不可导点的方法.令0)(*=∇x f ,求出平稳点.如果)(*2x f ∇是正定的,则*x 是UMP 的严格局部最优解.若)(x f 在n R 上是凸函数,则是整体最优解.在求解0)(*=∇x f 这n 维方程组比较困难时，就用最优化算法——迭代法.本节将介绍最速下降法,牛顿法,共轭方向法,坐标轮换法,变尺度法.这些算法就是用不同的方法来选择搜索方向k p 而得到的.当然kp 必须是下降方向.定理7.7 设R R f n→:，在点x 处可微,若存在nR p ∈,使0)(<∇p x f T,则向量p是f 在x 处的下降方向.证明:)(x f 可微(在x 处)，由泰勒展开式,有 ||)(||)()()(p p x f x f p x f Tλολλ+∇+=+ ,0,0)(><∇λp x f T0)(<∇∴p x f Tλ），（当δλδ0∈∃∴时,有0||)(||)(<+∇p p x f Tλολ),0()()(δλλ∈∀<+∴x f p x fp ∴是)(x f 在点x 的下降方向. 证毕.一、最速下降法最速下降法又称梯度法,选择负梯度方向作为目标函数值下降的方向,是比较古老的一种算法，其它的方法是它的变形或受它的启发而得到的，因此它是最优化方法的基础. 基本思想:迭代法求解无约束最优化(5.9)问题的关键是求下降方向kp .显然最容易想到的是使目标函数值下降速度最快的方向.从当前点kx 出发,什么方向是使)(x f 下降速度最快呢? 由泰勒展开知:||)(||)()()(k k T k k k k p p x f p x f x f λολλ+∇-=+-略去λ的高阶无穷小项,取)(kkx f p -∇=时,函数值下降最多.而)(kx f ∇为)(x f 在kx 处的梯度，所以下降方向kp 取为负梯度方向时,目标函数值下降最快.最速下降法:①. 取初始点0x ,允许误差0>ε,令0:=k ②. 计算)(kkx f p -∇=③. 若ε<||||k p ,停止，点k x 为近似最优解.否则进入④.④. 求 k λ,使)(min )(0kk k k k p x f p x f λλλ+=+≥ ⑤. 令kk k k p x xλ+=+1，1:+=k k ，返回②例7.4 用最速下降法求解下列无约束优化问题1222121225),(m in x x x x x f -+=取初始点Tx )2,2(0= 终止误差 610-=ε解：很显然，该问题的整体最优解为Tx )0,1(*=⎪⎪⎭⎫⎝⎛-=∇215022)(x x x f ，令0,10)(21==⇒=∇x x x f易验证)(*2x f ∇是正定的， ∴是整体最优解. 下面用最速下降法求解T T x x x f x f x f )50,22(),()(2121-=∂∂∂∂=∇ T x )2,2(0=T x f )100,2()(0=∇∴取Tp )100,2(0-=由⎪⎪⎭⎫⎝⎛--=⎪⎪⎭⎫ ⎝⎛--+⎪⎪⎭⎫ ⎝⎛=+λλλλ10022210022200p x4)22(2)1002(25)22()(2200+---+-=+λλλλp x f得0)1002(5000)22(4=----=λλλd df020007679.0500008100080==⇒λ⎪⎪⎭⎫⎝⎛-=⎪⎪⎭⎫ ⎝⎛--+⎪⎪⎭⎫ ⎝⎛=+=0007679.0959984642.11002020007679.0220001p x x λ重复上述过程得 Tx )01824717.0,009122542.1(2=789850288.0)(,078282.0)(,100)(21-=-==x f x f x f图7-1从图7-1可知，{}kx 随着迭代次数的增加，越来越接近与最优解)0,1(，同时也看到，随着迭代次数的增加，收敛速度越来越慢.极小点附近沿着一种锯齿形前进，即产生“拉锯”现象：{}kx沿相互正交的方向小步拐进，趋于最优解的过程非常缓慢.这种现象怎样解释？如何克服？在求k λ时，由于)()(kkp x f λλϕ+=，求导得0)('=λϕ，k λ是)(λϕ的极小点.故有0)()('=⋅+∇=k T k k k k p p x f λλϕ，即0)(=⋅+∇kk k k p p x f λ，又)(11++-∇=k k x f p，即0)(1=⋅+k T k p p 或0)()(1=∇⋅∇+k T k x f x f .也就是最速下降法相邻两个搜索方向是彼此正交的.因此最速下降法是局部下降最快，但不一定是整体下降最快的.在实际应用中一般开始几步用最速下降法，后来用下面介绍的牛顿法.这样两个算法结合起来，求解速度较快.二、牛顿法 基本思想：考虑无约束优化问题(5.9))(min x f nRx ∈ 由前面的讨论知，若能解出方程组0)(=∇x f ，求出平稳点*x ，则可验证*x 是否为极值点.由于0)(=∇x f 难以求解.如果)(x f 具有连续的二阶偏导数,则考虑用)(x f 在点*x 二阶泰勒展开式条件替代)(x f ∇，即由22||)(||))(()(21)()()()(k k k T k k T k k x x x x x f x x x x x f x f x f -+-∇-+-∇+=ο))(()(21)()()()()(2kk T k k T k k x x x f x x x x x f x f x g x f -∇-+-∇+=≈⇒令0))(()()()(2=-∇+∇=∇≈∇kk k x x x f x f x g x f)())((121k k k k x f x f x x ∇∇-=⇒-+即从kx 出发，搜索方向为)())((12kkkx f x f p ∇∇-=-，步长恒为1，得到下一个迭代点1+k x.牛顿法：①． 选取初始点0,0=：k x ，精度0>ε ②． 计算)(kx f ∇，如果ε≤∇||)(||kx f ，计算终止.否则计算)(2kx f ∇，求出搜索方向)())((12kk k x f x f p ∇∇-=-. ③． 令1:,1+=+=+k k p x x k k k ，返回②.例7.5 考虑无约束问题22122214)(m in x x x x x f -+=试分别取初始点（1）T x )1,1(0=，（2）T x )4,3(0=（3）Tx )0,2(0=，取精度要求310-=ε，用牛顿法求解.解：⎪⎪⎭⎫ ⎝⎛--=∇212211228)(x x x x x x f ⎪⎪⎭⎫⎝⎛---=∇22228)(1122x x x x f (1) 取Tx )1,1(0=有Tx f )1,6()(0=∇ ε>=∇0828.6||)(||0x f⎪⎪⎭⎫⎝⎛--=∇2226)(02x fT x f x f p )2500.2,7500.1()())((01020--=∇⋅∇-=-Tp x x )2500.1,7500.0(01--=+= 重复计算结果得表7-2.ε<=0||)(||4x f T x )0,0(4=∴为近似最优解.实际上，该问题最优解为**)0,0(=x(2) 取Tx )4,3(0=,同上计算,得TT x x x )4,8284.2(,)4,8333.2(),4,3(21===有ε<=∇=∇=∇0||)(||,1667.0||)(||,1||)(||210x f x f x f ，这一迭代结果收敛到)(x f 的鞍点T)4,22(.(3) 取Tx )0,2(0=T x f )4,16()(0-=∇ ⎪⎪⎭⎫⎝⎛--=∇2448)(02x f0)(02=∇x f , 即)(02x f ∇不可逆，所以无法求得点1x .牛顿法的主要缺点:（1） 该法的某次迭代反而使目标函数值增大(见上例).（2） 初始点0x 距极小点*x 较远时,产生的点列{}kx可能不收敛,还会出现)(*2x f ∇的奇异情况.（3） )(*2x f ∇的逆矩阵计算量大. 牛顿迭代法的主要优点:当目标函数)(x f 满足一定条件,初始点0x 充分接近极小点*x 时,由牛顿法产生的点列{}kx 不仅能够收敛到*x,而且收敛速度非常快.实际应用中常将最速下降法和牛顿法结合起来使用.牛顿法的改进:上面介绍的牛顿法中,kx 处的搜索方向为)())((12kkkx f x f p ∇∇-=-,步长恒为 1.若通过一维搜索来取最优步长k λ,可防止在迭代中步长恒为1时标目标函数值增大的可能. 改进的牛顿法:①. 取初始点nR x ∈0,允许误差0:,0=>k ε.②. 检验是否满足ε<∇||)(||kx f ,若是,迭代停止,得到k x 为近似最优解.否则进入③.③. 令)())((12kk k x f x f p ∇∇-=-.④. 求k λ,使)()(min kk k k k p x f p x f λλλ+=+. ⑤. 令k k k k p x x λ+=+1,令1+=k k ：转②.三、坐标轮换法既然求解非线性规划问题的迭代法是给出初始点0x ,求出一个方向kp ,根据kp 确定步长k λ,使k k k k p x xλ+=+1，如果1+k x 满足某精度要求则停止，否则继续找方向1+k p .显然构造出搜索方向有一定的困难，能否按既定的搜索方向寻找最优解,省去找搜索方向kp 呢?在最速下降法中我们看到相邻两个搜索方向kp 和1+k p是正交的.我们知道在n 维欧氏空间中坐标轴向量n εεε,,,21 是正交的，可否选坐标轴向量为搜索方向kp 为呢?回答是肯定的，这样我们得到了坐标轮换法.基本思想:从1x 出发,取11ε=p ,沿1p 进行一维搜索得到1112p x x λ+=.若2x 满足精度要求,则停止.否则取22ε=p ,2223p x x λ+=.如此继续,, 取n n n n n n p x x p λε+==+1,，若1+n x 满足精度要求,则停止.否则令11ε=+n p ,1112+++++=n n n n p x x λ,如此反复连续,可以求出近似最优解.坐标轮换法的缺点是收敛速度较慢,搜索效率较低,但基本思想简单,沿坐标轴的方向进行搜索.四、共轭方向法和共轭梯度法共轭方向法是一类方法的总称,它原来是为求解目标函数为二次函数的问题而设计的.这类方法的特点是:方法中的搜索方向是与二次函数的系数矩阵Q 有关的所谓共轭方向，用这类方法求解n 元二次函数的极小化问题最多进行n 次一维搜索便可以得到极小点.由于可微的非二次函数在极小点附近的性态近似于二次函数,因此这类方法也用于求可微的非二次函数的UMP 问题.定义7.14 设Q 为n n ⨯对称正定矩阵,如果0=Qy x T称n 维向量x 和y 关于Q 共轭.定义7.15 设k p p p ,,,21 为nR 中一组向量, Q 是一个n n ⨯对称正定矩阵.如果k j i j i Qp p Qp p i T i j T i ,,2,1,,,0,0 =≠≠=,称k p p p ,,,21 为Q 共轭向量组,也称它们为一组Q 共轭方向.当E Q =(单位矩阵)时,为正交方向.定理7.8 若k p p p ,,,21 为矩阵Q 共轭向量组,则它们必线性无关. 证明: 若存在k l l l ,,,21 ,使011=++k k p l p l ,则对任一k j ,,2,1 =,由 0)(11===∑∑==j T j j ki j T j iki iiT jQp p l Qp pl p l Q p又0>j Tj Qp p , 0=∴j l∴ k p p p ,,,21 线性无关. 证毕.由高等代数知识可知, Q 共轭向量组中最多包含n 个向量, n 是向量的维数.反之,可以证明,由n 维空间的任一组基出发可以构造出一组Q 共轭方向11,,,-n pp p .前面我们已经讲述了坐标轮换法，在n 维欧几里德空间中, n εεε,,,21 是一组线性无关的正交向量.从0x 出发,依次使用n εεε,,,21 作为下降方向进行一维精确搜索来确定n x x x ,,,21 ,重复进行得点列{}k x ,最终求得满足精度要求的最优解.刚才我们引进了共轭向量组11,,,-n p p p ,又证明了它们的线性无关性,那么是否可以用这共轭向量组像坐标轮换法一样,从0x 出发依次用11,,,-n pp p 作为下降方向来确定{}kx,最终求得近似最优解?回答是肯定的.这个方法称为共轭方向法.共轭方向法适合任何可微凸函数,且对于二次函数极值)(min x f x p Qx x T T+=21特 别有效.下面的定理告诉我们用共轭方向法求解二次函数的极值，经过n 次迭代就能求得最优解.定理7.9 设Q 为n n ⨯对称正定矩阵,函数x p Qx x x f T T+=21)(,又设 110,,,-n p p p 为一组Q 共轭向量组,且每个i p 是(下面形成的)i x 点处的下降方向.则由任一点0x 出发,按如下迭代至多n 步后收敛,k k k k p x xλ+=+1,这里k λ满足)(m i n )(0k k k k k p x f p x f λλλ+=+>.证明: 欲证至多n 步收敛,即证0)(=∇nx f .下证)(nx f ∇和11,,,-n pp p 正交.p Qx x f +=∇)( p Qx x f kk+=∇∴)( p p x Q p Qx xf k k k k k ++=+=∇++)()(11λkk k k k k Qp x f p Qp Qx λλ+∇=++=)( =+∇=∇---111)()(n n n n Qpx f x f λ 11111)(--++++++∇=n n k k k Qp Qp xf λλQ p Q p x f x f Tn n T k k T k T n )()()()(11111--++++++∇=∇λλkT n n k T k k k T k k T n Qp p Qp p p x f p x f )()()()(11111--++++++∇=∇λλ000+++= )2,,2,1,0(-=n k 又0)(1=∇-n Tn px f0)(=∇∴kT n p x f )1,,1,0(-=n k)(nx f ∇∴和11,,,-n pp p 正交， 又110,,,-n pp p 线性无关.0)(=∇∴nx fnx ∴是问题的最优解. 证毕. 共轭方向法具有二次有限终止性. 由于共轭方向组11,,,-n p p p 的取法有很大的随意性,用不同方式产生一组共轭方向就得到不同的共轭方向法.如果利用迭代点处的负梯度向量为基础产生一组共轭方向,这样的方法叫共轭梯度法.下面对二次函数讨论形成Q 共轭梯度方向的一般方法,然后引到求解无约束化问题上.任取初始点n R x ∈0,若0)(0≠∇x f ,取)(0x f p -∇=,从0x 点沿方向0p 进行一维搜索 ,求得0λ.令0001p x x λ+=,若0)(1=∇x f ,则已获得最优解1*x x =.否则，取0011)(p x f p υ+-∇=,其中0υ的选择要使1p 和0p 关于Q 共轭，由0)(01=Qp p T ,得0100)()()(Qp p x f Q p T T ∇=υ一般地,若已获得Q 共轭方向kp p p ,,,1和依次沿它们进行一维搜索的得到的点列110,,,+k x x x ,若0)(1=∇+k x f ,则最优解为1*+=k x x ,否则∑=+++-∇=ki i i k k p xf p011)(α为使1+k p 和11,,,-k pp p 是共轭,可证0110====-k ααα所以有 k k k k p x f pυ+-∇=++)(11又1+k p和kp 是Q 共轭的.有0)(1=+k Tk Qp p，得kT k k T k k Qpp x f Q p )()()(1+∇=υ 2,,2,1,0-=n k 进一步可得k υ221||)(||||)(||k k x f x f ∇∇=+ 2,,1,0-=n k综合起来得 Fletcher-Reeves 公式2)21110||(||||)(||)()(k k k k k k k x f x f p x f px f p ∇∇=+-∇=-∇=+++υυ 2,,2,1,0-=n k （7.10）共轭梯度法： ①． 选取初始点0x ，给定终止误差0>ε. ②． 计算)(0x f ∇，若ε≤∇||)(||0x f ,停止迭代，输出0x .否则进行③.③． 取)(0x f p -∇=，令0:=k④． 求k λ，)(min )(0kkkk kp x f p x f λλλ+=+≥，令k k k k p x xλ+=+1⑤． 计算)(1+∇k xf ，若ε≤∇+||)(||1k x f ，停止迭代，1*+=k x x 为最优解.否则转⑥.⑥． 若n k =+1，令nx x =:0，转③(已经完成一组共轭方向的迭代,进入下一轮)否则转⑦ ⑦． 取kk k k p xf pυ+-∇=++)(11,其中221||)(||||)(||k k k x f x f ∇∇=+υ,令1:+=k k ,转④当)(x f 是二次函数时上述共轭梯度法至多进行n 步可求得最优解.当)(x f 不是二次函数,共轭梯度法也可以构造11,,,-n p p p ,但已不具有有限步收敛的性质，于是和坐标轮换法一样经过一轮迭代后,采用重新开始的技巧.上述共轭梯度法就是带有再开始技巧的F-R 法.例7.6 用F-R 法求下面问题 2212121252),(m in x x x x x f +-=取初始点T x )2,2(0=,终止误差为610-=ε解:在例7.4中已得Tx f p )100,2()(0-=-∇= Tx )0007679.0,959984642.1(1-= Tx f )038395.0,919969284.1()(1-=∇000368628.010004687756228.3||)(||||)(||20210==∇∇=x f x f υ ⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎭⎫ ⎝⎛--+⎪⎪⎭⎫ ⎝⎛-=+-∇=0015322.092070654.11002000368628.0038395.0919969284.1)(0011p x f p υ⎪⎪⎭⎫ ⎝⎛+--=+0015322.00007679.092070654.1959984642.111λλλp x0378228399.7687703443.3)(11=+-=+λλλd p x df499808794.01=∴λ⎪⎪⎭⎫ ⎝⎛≈⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎭⎫⎝⎛⨯+--⨯+=+=010********.0999998622.00015322.0499808794.00007679.0)92070654.1(499808794.0959984642.11112p x x λε<=∇0||)(||2x f , ∴最优解⎪⎪⎭⎫⎝⎛==012*x x .五、变尺度法当初始点为)(x f 的其极值点附近时牛顿法收敛速度很快，但缺点是需计算)(2kx f ∇的逆矩阵，在实际问题中目标函数往往相当复杂，计算二阶导数的工作量或者太大或者不可能，在x 的维数很高时，计算逆矩阵也相当费事.如果能设法构造另一个矩阵kH ，用它来逼近二阶导数矩阵的逆矩阵12))((-∇kx f 则可避免上述问题.下面就来研究如何构造12))((-∇kx f 的近似矩阵kH .我们希望：每一步都能以现有的信息来确定下一个搜索方向，每做一次迭代，目标函数值均有所下降，这些近似矩阵最后应收敛于最优解处的海赛矩阵的逆矩阵12))((-∇kx f .p Qx x f xp Qx x x f T T+=∇+=)(21)(考虑于是 )]()([)()()(11111k k k k k k k k x f x f Q x x x x Q x f xf ∇-∇=-⇒-=∇-∇+-+++当)(x f 是非二次函数时，令)]()([111k k k k k x f x f H x x ∇-∇=-+++ （7.11）称为拟牛顿条件.显然，我们构造出来的近似矩阵k H 必须满足上述拟牛顿条件及递推性：k k k H H H ∆+=+1，这里k H ∆称为矫正矩阵.记 k k k kk k x x x x f x f G -=∆∇-∇=∆++11)()( 有 kk k k k k G H H G H x ∆∆+=∆=∆+)(1 .变尺度法即DEP 法是由Davidon 首先提出，后来又被Fletcher 和Powell 改进的算法.记kk T k kT k k k k T k T k k k k kk T k kT k k k k T k T k k kG H G HG G H x G x x H H G H G H G G H x G x x H ∆∆∆∆-∆∆∆∆+=∆∆∆∆-∆∆∆∆=∆+)()()()()()()()(1 （7.12）容易验证1+k H 满足拟牛顿条件，称上式为DEP 公式.变尺度方法计算步骤：（１） 给出初始点nR x ∈0，允许误差0>ε.（２） 若ε<∇||)(||0x f ，停止，0x 为近似最优解；否则转下一步.（３） 取I H =0（单位矩阵），0=：k . （４） )(kk k x f H p ∇-=（５） 求k λ，使)(min )(0kk k k k p x f p x f λλλ+=+≥. （６） 令kk k k p x xλ+=+1（７） 若ε<∇+||)(||1k xf ，1+k x 为最优解，停止；否则当1-=n k 时，令n x x =：0转（３）.（即迭代一轮n 次仍没求得最优解，以新的0x 为起点重新开始一轮新的迭代）.k k k k k kx x x x f xf G n k -=∆∇-∇=∆-<++11),()(1时，令当.计算kk T k kT k k k k T k T k k kk G H G H G G H x G x x H H∆∆∆∆-∆∆∆∆+=+)()()()(1，令1+=k k ：，转（４）. §4 约束极值的最优化方法考虑(MP)问题：0)(0)(..)(min =≥x H x g t s x f （7.13）其中Tl T m x h x h x h x g x g x g ))(,),(()(,))(,),(()(11 ==是向量函数.显然 0)(0)(0)(≥-≥⇔=x h x h x h ， 于是（MP ）问题可以写为：Tm x g x g x g x g t s x f ))(,),(()(0)(..)(min 1 =≥ （7.14）一、积极约束设0x 是（MP ）问题（5.14）的一个可行解.对0)(0≥x g i 来说，在点0x 有两种情况：或者0)(0>x g i ，或者0)(0=x g i .0)(0>x g i 时，则0x 不在0)(0=x g i 形成的边界上，称这一约束为0x 的非积极约束；0)(0=x g i 时，0x 处于由0)(0≥x g i 这个约束条件形成的可行域边界上，当0x 有变化时就不满足0)(0=x g i 的条件，所以称为积极约束，记为：{}()|()0,1i I x i g x i m ==≤≤.定义7.16 设x 满足约束条件0)(0≥x g i ),,1(m i =，0)(|{)(==x g i x I i ，}m i ≤≤1，如果)(x g i ∇，)(x I i ∈线性无关，则称点x 是约束条件的一个正则点.二、可行方向、下降方向的代数条件前面我们已经给出可行方向和下降方向的定义，下面给出其代数条件.可行方向：设K 是（MP ）问题（5.14）的可行域，K x ∈,0,≠∈p R p n.若存在00>λ使得],0[0λλ∈时有K p x ∈+λ，称p 为x 点处的一个可行方向.由泰勒公式：||)(||)()()(p p x g x g p x g T i i i λολλ+∇+=+当x 为)(x g i 的积极约束时，有0)(=x g i .只要0>λ足够小，)(p x g i λ+和p x g T i )(∇λ同号，于是当0)(>∇p x g T i 时有0)(≥+p x g i λ )(x I i ∈.当x 为)(x g i 的非积极约束时，有0)(>x g i .由)(x g i 的连续性，当0>λ足够小时，由保号性知 0)(≥+p x g i λ )(x I i ∉.所以只要 0)(>∇p x g T i ，)(x I i ∈就可保证0)(≥+p x g i λ，于是p 为x 点处的一个可行方向.称0)(>∇p x g T i ，)(x I i ∈ 为p 在点x 处是可行方向的代数条件.下降方向：设K 是（MP ）问题的可行域，K x ∈,0,≠∈p R p n.若存在00>λ使得],0[0λλ∈时，有)()(x f p x f <+λ，称p 为x 点处的一个下降方向.由泰勒公式：||)(||)()()(p p x f x f p x f Tλολλ+∇+=+.当λ足够小时，只要0)(<∇p x f T，有)()(x f p x f <+λ. 称0)(<∇p x f T为p 在x 点处的一个下降方向的代数条件.三、可行下降方向设K 为（MP ）问题（5.14）的可行域，K x ∈，若存在0,≠∈p R p n，p 既是x 点处的下降方向又是可行方向，则称p 为点x 处的可行下降方向.定理7.10 考虑非线性规划问题（5.14），K x ∈，),,1)()(m i x g x f i =（和在x点处可微，若*x 是局部极小点，则x 点处必不存在可行下降方向，即不存在p 同时满足：⎪⎩⎪⎨⎧∈>∇<∇)(0)(0)(x I i p x g p x f Ti T证明：反证法，设局部极小点x 处存在可行下降方向p ，于是1λ∃，当],0[1λλ∈时有)()(x f p x f <+λ，又p 为可行方向．2λ∃∴当],0[2λλ∈时K p x ∈+λ，这与x 是。
数值分析实验报告_清华大学_非线性方程的解法
非线性方程的解法实验1.算法设计与比较问题提出：非线性方程组的求解方法很多，基本的思想是线性化。
不同的方法效果如何，要靠计算的实践来分析、比较。
实验内容：考虑算法(1)牛顿法(2)拟牛顿法分别编写它们的matlab程序。
实验要求：(1)用上述方法，分别计算两个例子。
在达到精度相同的前提下，比较迭代次数、浮点运算次数和CPU时间等。
1.1程序清单为使用flops统计浮点运算次数，使用MATLAB5.3版本%f1.m原函数f1function y=f(x)y(1)=12*x(1)-x(2)^2-4*x(3)-7;y(2)=x(1)^2+10*x(2)-x(3)-8;y(3)=x(2)^3+10*x(3)-8;end%ff1.m原函数f1的雅克比矩阵function y=ff(x)y(1,:)=[12,-2*x(2),-4];y(2,:)=[2*x(1),10,-1];y(3,:)=[0,3*x(2)^2,10];end%f1.m原函数f2function y=f2(x)y(1)=3*x(1)-cos(x(2)*x(3)) -1/2;y(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;y(3)=exp(-x(1)*x(2))+20*x(3)+1/3*(10*pi-3);end%ff2.m原函数f2的雅克比矩阵function y=ff2(x)y(1,:)=[3,x(3)*sin(x(2)*x(3)),x(2)*sin(x(2)*x(3))];y(2,:)=[2*x(1),-2*81*(x(2)+0.1),cos(x(3))];y(3,:)=[-x(2)*exp(-x(1)*x(2)),-x(1)*exp(-x(1)*x(2)),20]; end%牛顿法(以第一个方程组为例)clear;x0=[0,0,0]';n=10;tol=1e-6;x(:,1)=x0;i=1;u=[1,1,1]';tic;while (norm(u)>tol*norm(x(:,i))&(i<n))A=ff1(x(:,i));b=f1(x(:,i))';u=-A\b;x(:,i+1)=x(:,i)+u;i=i+1;end;x(:,i)iter=i-1t=toc%拟牛顿法(以第一个方程组为例)clear;x0=[0,0,0]';n=10;tol=1e-6;x(:,1)=x0;i=1;p=[1,1,1]';A=ff1(x(:,1));tic;while (norm(p)>tol*norm(x(:,i))&(i<n))x(:,i+1)=x(:,i)-A\f1(x(:,i))';p=x(:,i+1)-x(:,i);q=f1(x(:,i+1))'-f1(x(:,i))';A=A+(q-A*p)*p'/norm(p,2)^2;i=i+1;end;iter=i-1t=tocx(:,i)1.2运行结果1.2.1第一个方程组精确解为*T =(0.886020214719037, 0.796444775323146, 0.749479574122230)x 取最大迭代次数n=5000，相对误差限Tol=1e-6 (1)取()(0)1,1,1x T=牛顿迭代法迭代3次收敛，浮点运算次数为440，每次迭代平均浮点运算次数为147，CPU 耗时t =0(s)拟牛顿法迭代4次收敛，浮点运算次数为1048，每次迭代平均浮点运算次数为262，CPU 耗时t =0(s)(2)取()(0)000x T =，， 牛顿迭代法迭代4次收敛，浮点运算次数为510，每次迭代平均浮点运算次数为128，CPU 耗时t =1.600e-002(s)拟牛顿法迭代6次收敛，浮点运算次数为1493，每次迭代平均浮点运算次数为248，CPU 耗时t =1.50e-002(s)(3)取()(0)50,5050x T=，牛顿迭代法迭代15次收敛，浮点运算次数为2118，每次迭代平均浮点运算次数为141，CPU 耗时t =1.600e-002(s)拟牛顿法迭代338次收敛，浮点运算次数为88454，每次迭代平均浮点运算次数为262，CPU 耗时t =3.100e-002(s)1.2.2第二个方程组精确解为*T =(0.886020214719037, 0.796444775323146, 0.749479574122230)x 取最大迭代次数n=5000，相对误差限Tol=1e-6(1)取()(0)000x T=，， 牛顿迭代法迭代5次收敛，浮点运算次数为776，每次迭代平均浮点运算次数为155.2，CPU 耗时t =0(s)拟牛顿法迭代6次收敛，浮点运算次数为1635，每次迭代平均浮点运算次数为273，CPU 耗时t =0(s)(2)取()(0)888x T=，， 牛顿迭代法迭代9次收敛，浮点运算次数为1519，每次迭代平均浮点运算次数为169，CPU 耗时t =0(s)拟牛顿法迭代21次收敛，浮点运算次数为5924，每次迭代平均浮点运算次数为282，CPU 耗时t =1.600e-002(s)(3)对于离精确解更远的初值（如()(0)101010x T=，，），在计算中会出现奇异或接近奇异的矩阵，计算结果误差很大或计算根本无法进行下去。
实验九 非线性规划
min f ( X ) h i ( X ) 0 , i 1, 2 , , m g ( X ) 0 , j 1, 2 , , l j
@for(num_j(j):gj(x)>=0;); !不 等 式 约 束 条 件 @for(num_k(k):x(k)>=0;); END
上次课的内容回顾：
(1)带有约束非线性规划的一般模型:
min f ( X ) h i ( X ) 0 , i 1, 2 , , m g ( X ) 0 , j 1, 2 , , l j
(2)介绍了非线性规划的最优性条件； (3)非线性规划的求解方法和步骤； (4)非线性规划的LINGO模型。
[OBJ]max=20*x(1)+16*x(2) -2*x(1)^2-x(2)^2-(x(1)+x(2))^2; x(1)+x(2)<=5; @for(num_i(i):x(i)>=0;); END
3
2012年11月7日
问题的LINGO模型： MODEL: sets: num_i/1,2/:x; Endsets [OBJ]max=20*x(1)+16*x(2) -2*x(1)^2-x(2)^2-(x(1)+x(2))^2; x(1)+x(2)<=5; @for(num_i(i):x(i)>=0;); END
1
2012年11月7日
上次课的内容回顾：
MODEL: sets: num_i/1..m/; num_j/1..L/; num_k/1..n/:x0,x; endsets init: x0=x0(1),x0(2),… ,x0(n); endinit [OBJ]min=f(x); !目 标 函 数 的 表 达 式 !等 式 约 束 条 件 @for(num_i(i):hi(x)=0;);
第九章非线性规划
资料范本本资料为word版本，可以直接编辑和打印，感谢您的下载第九章非线性规划地点：__________________时间：__________________说明：本资料适用于约定双方经过谈判，协商而共同承认，共同遵守的责任与义务，仅供参考，文档可直接下载或修改，不需要的部分可直接删除，使用时请详细阅读内容1. 非线性规划我们讨论过线性规划,其目标函数和约束条件都是自变量的线性函数。
如果目标函数是非线性函数或至少有一个约束条件是非线性等式(不等式),则这一类数学规划就称为非线性规划。
在科学管理和其他领域中，很多实际问题可以归结为线性规划，但还有另一些问题属于非线性规划。
由于非线性规划含有深刻的背景和丰富的内容，已发展为运筹学的重要分支，并且在最优设计,管理科学,风险管理,系统控制,求解均衡模型,以及数据拟合等领域得到越来越广泛的应用。
非线性规划的研究始于三十年代末,是由W.卡鲁什首次进行的，40年代后期进入系统研究，1951年H.W.库恩和A.W.塔克提出带约束条件非线性规划最优化的判别条件,从而奠定了非线性规划的理论基础,后来在理论研究和实用算法方面都有很大的发展。
非线性规划求解方法可分为无约束问题和带约束问题来讨论，前者实际上就是多元函数的极值问题，是后一问题的基础。
无约束问题的求解方法有最陡下降法、共轭梯度法、变尺度法和鲍威尔直接法等。
关于带约束非线性规划的情况比较复杂，因为在迭代过程中除了要使目标函数下降外，还要考虑近似解的可行性。
总的原则是设法将约束问题化为无约束问题；把非线性问题化为线性问题从而使复杂问题简单化。
求解方法有可行方向法、约束集法、制约函数法、简约梯度法、约束变尺度法、二次规划法等。
虽然这些方法都有较好的效果，但是尚未找到可以用于解决所有非线性规划的统一算法。
1.1 非线性规划举例[库存管理问题] 考虑首都名酒专卖商店关于啤酒库存的年管理策略。
假设该商店啤酒的年销售量为箱,每箱啤酒的平均库存成本为元,每次订货成本都为元。
【清华】Exp-09
min f (x) = 1 xT Hx + cx 2
s.t. Ax ≤ b
(1)
min f (x) = 1 xT Hx+ cx
2
(2)
s.t. a x = b ,
j
j
j∈J 1
(a j是 A的第 j列)
基本 原理
• 若x为(1)的最优解，则它也是(2)的最优解
• 若x为(1)的可行解，又是(2)的KKT点，且 L—乘子非负，则它必是(1)的最优解。
1
2
s.t. g (x) = x2 + x2 −10≤ 0
1
1
2
g (x) = x + x −4 ≤ 0
2
1
2
x2
g (x) = x ≥ 0
3
2
∇f
Q
(7,3) • 最优解在P(3,1)取得
-∇g2 -∇g1
-∇g1
P
λ1 =1, λ2 = 2, λ3 = 0, ∇f (P) +∇g1(P) +2∇g2(P) = 0
QP子 问题
min 1 d T G d + ∇f (x )T d
2
k
k
s.t.
∇h (x )T d + h (x ) = 0,
ik
ik
i = 1,Lm
∇g (x )T d + g (x ) ≤ 0,
jk
jk
j = 1,Ll
x 是第k 次迭代的初始点，G 是海赛阵∇2 L 的近似。
k
k
将最优解dk作为迭代的搜索方向，令 xk+1 = xk +αkdk
λ1,L λl
≥0
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验9 非线性规划实验目的：1）掌握用matlab优化工具箱解非线性规划的方法2）练习建立实际问题的非线性规划模型实验内容：4.某公司将3种不同含硫量的液体原料（分别记为甲、乙、丙）混合生产两种产品（分别记为A，B）.按照生产工艺的要求，原料甲、乙必须首先倒入混合池中混合，混合后的液体再分别于原料丙生产A，B.已知原料甲、乙、丙的含硫量分别是3%，1%，2%，进货价格分别为6千元/t，16千元/t，10千元/t；产品A,B的含硫量分别不能超过2.5%，1.5%，售价分别为9千元/t，15千元/t.根据市场信息，原料甲、乙、丙的供应量都不能超过500t;产品A，B的最大市场需求量分别为100t，200t.（1）应如何安排生产？（2）如果产品A的最大市场需求量增长为600t，应如何安排生产？（3）如果乙的进货价格下降为13千元/t，应如何安排生产？分别对（1）、（2）两种情况进行讨论.解：（1）问题的建模设利用x1吨甲，x2吨乙，x3吨丙制造y1吨A；利用x2吨甲，x4吨乙，x6吨丙制造y2吨B；总收益是z千元。
则有以下方程与不等式：质量守恒：y1=x1+x3+x5y2=x2+x4+x6总收益：z=9y1+15y2-6(x1+x2)-16(x3+x4)-10(x5+x6)化简得：z=3x1+9x2+3x3+9x4-x5+5x6含硫量约束：3%x1+1%x3+2%x5≤2.5%y13%x2+1%x4+2%x6≤1.5%y2化简得：0.5 x1-1.5x3-0.5x5≤01.5x2-0.5x4+0.5x6≤0供应量约束：（x1+x2）,(x3+x4),(x5+x6)≤500需求量约束：y1≤100；y2≤200化简得：x1+x3+ x5≤100x2+x4+ x6≤200甲乙混合，比例相同：x1 3=x24整理得：x1x4-x2x3=0;模型的求解：该问题是一个带约束非线性规划问题，编写源程序如下：M文件：函数文件：function z = lab94fun( x)z=-(3*x(1)+9*x(2)-7*x(3)-x(4)-x(5)+5*x(6));end非线性约束条件文件：function [ c,ceq ] = lab94con( x )c=0;ceq=x(1)*x(4)-x(2)*x(3);end主文件：A=[0.5 0 -1.5 0 -0.5 00 1.5 0 -0.5 0 0.51 1 0 0 0 00 0 1 1 0 00 0 0 0 1 11 0 1 0 1 00 1 0 1 0 1 ]b=[0 0 500 500 500 100 200]'x0=[20 20 20 20 30 30]; %已验证在可行域中v1=[0 0 0 0 0 0][x,z,ef,out,lag,grad,hess]=fmincon(@lab94fun,x0,A,b,[],[],v1,[],@lab94con); xz运算结果为：x =0.0000 0 -0.0000 100.0000 0 100.0000-400因此，此时应购买100吨乙，100吨丙来生产200吨B，总共收益是400千元（2）问题的建模：修改x1+x3+ x5≤100为：x1+x3+ x5≤600，其余不变。
模型的求解：主文件：将b变为：b=[0 0 500 500 500 600 200]'在实际运行时发现，在初值为x0=[150 10 10 10 150 10]（也有其他初值）时，总收益最大。
结果为：x =300.0000 0.0000 0 -0.0000 300.0000 0.0000z =-600因此，这时应该购入100吨甲和100吨丙来生产A，总收益是600千元。
（3）问题的建模：总收益变为：z=3x1+9x2-4x3+2x4-x5+5x6其余不变。
模型的求解：将函数文件修改为：z=-(3*x(1)+9*x(2)-4*x(3)+2*x(4)-x(5)+5*x(6));其余不变。
对（1），结果为：（修改初值是x0=[30 10 10 100 30 10]）x =0.0000 50.0000 -0.0000 150.0000 0 0z =-750.0000对（2），结果为：（修改初值是x0=[30 10 10 100 30 10]）0.0000 50.0000 -0.0000 150.0000 0 0z =-750.0000可见，二者结果相同。
因此，应购买50吨甲，150吨乙，生产200吨B，总收益是750千元。
（本题初值对结果有一定影响，有时初值不同、结果不同，因此应该多选初值进行尝试。
）9.8.美国某三种股票（A，B，C）12年（1943-1954年）的价格（已经包括了分红在内）每年的增长情况如下表所示（表中还给出了相应年份的500种股票的价格指数的增长情况）.例如，表中第一个数据1.300倍，即收益为30%，其余数据的含义以此类推.假设你在1955年时有一笔资金准备投资这三种股票，并期望年收益率至少达到15%，那么你应当如何投资？此外，考虑一下问题：（1）当期望的年收益率在10%~100%变化时，投资组合和相应的风险如何变化？（2）假设除了上述三种股票外，投资人还有一种无风险的投资方式，如购买国库券.假设国库券的年收益率为5%，如何考虑该投资问题？（3）假设你手上目前握有的股票比例为：股票A占50%，B占35%，C占15%.这个比例与你得到的最优解可能有所不同，但实际股票市场上每次股票买卖通常总有交易费，例如按交易额的1%收取交易费，这时你是否仍需要对手上的股票进行买卖（换手），以便满足“最优解”的要求？年份股票A 股票B 股票C 股票指数1943 1.300 1.225 1.149 1.258997 1944 1.103 1.290 1.260 1.197526 1945 1.216 1.216 1.419 1.364361 1946 0.954 0.728 0.922 0.919287 1947 0.929 1.144 1.169 1.057080 1948 1.056 1.107 0.965 1.055012 1949 1.038 1.321 1.133 1.187925 1950 1.089 1.305 1.732 1.317130 1951 1.090 1.195 1.021 1.240164 1952 1.083 1.390 1.131 1.183675 1953 1.035 0.928 1.006 0.9901081954 1.176 1.715 1.908 1.526236 解：问题的建模：为方便建模及程序编写，先将各值计算出来。
令EA,EB,EC为各股票的期望，DA,DB,DC 为各股票的方差，rAB、rAC、rBC分别为每两支股票之间的相关系数。
可以计算得：EA=1.0891;EB=1.2137;EC=1.2346DA=0.0108;DB=0.0584;DC=0.0942rAB=0.4939;rAC=0.4097;rBC=0.7472同时可计算出每两种股票的协方差为：cov A,B=rAB DA DB=0.01240387cov A,C=rAC DA DC=0.01306782cov B,C=rBC DB DC=0.05542028令投资A,B,C三种股票每种所占的比例分别为x1、x2、x3因此，投资的总期望收益是：Z1=x1EA+x2EB+x3EC=1.0891x1+1.2137x2+1.2346x3投资总收益的方差为Z2=x12DA+x22DB+x32DC+2x1x2cov A,B+2x1x3cov A,C+2x2x3cov B,C=0.0108x12+0.0584x22+0.0942x32+0.01240387x1x2+0.01306782x1x3+0.05532028x2x3年收益率至少要达到15%：1.0891x1+1.2137x2+1.2346x3≥1.15其他约束：x1+x2+x3=1x1,x2,x3≥0模型的求解：源程序如下：函数文件：function f=lab98fun(x)f=0.0108*x(1)^2+0.0584*x(2)^2+0.0942*x(3)^2+0.01240387*x(1)*x(2)+0.01306782*x(1)*x(3)+0.05 532028*x(2)*x(3);end约束条件文件：function [c,ceq]=lab98con(x)c=1.15-1.0891*x(1)-1.2137*x(2)-1.2346*x(3);ceq=x(1)+x(2)+x(3)-1;end主程序：x0=[0.1 0.1 0.1];v1=[0 0 0];v2=[1 1 1];[x,fv,ef,out,lag,grad,hess]=fmincon(@lab98fun,x0,[],[],[],[],v1,v2,@lab98con);运行结果为：x =0.5394 0.2928 0.1678fv =0.0167所以选择投资A54%，B29%，C17%左右，可使方差最小，收益率的方差为1.67%。
由表中数据所示，当期望年收益率低的时候，多投资A股一些可以减小风险。
随着期望年收益率的不断增高，投资A股的份额不断减小，而B股和C股的份额在增大，同时风险也在增大，当期望年收益率更高时，B股的股份也下降，转而都投资C股，同时风险进一步加大。
由此可以看出，高收益总是伴随着高风险。
（2）问题的建模：只是要改变相应的约束条件。
ED=1.05 DD=0投资的总期望收益变为Z1=x1EA+x2EB+x3EC=1.0891x1+1.2137x2+1.2346x3+1.05x4投资总收益的方差不变其余约束条件变为x1+x2+x3+x4=1x1,x2,x3,x4≥0根据问题要求，变为1.0891x1+1.2137x2+1.2346x3+1.05x4≥1.15模型的求解：解为：x =0.2463 0.3268 0.1997 0.2272fv =0.0159因此，应投资的比例分别为A24.6%，B32.7%，C20%，D22.7%左右，收益率的方差为1.59%（3）问题的建模：需要改变收益的方程为Z1=x1EA+x2EB+x3EC=1.0891x1+1.2137x2+1.2346x3−0.01（0.5−x1+0.35−x2+0.15−x3）模型的求解：将收益的语句改为c=1.15-1.0891*x(1)-1.2137*x(2)-1.2346*x(3)+0.01*(abs(0.5-x(1))+abs(0.35-x(2))+abs(0.15-x(3)));运行后可得x =0.5314 0.3187 0.1500fv =0.0169因此，应投资的比例分别为A53.14%，B31.87%，C15%，收益率的方差为1.69%也就是适当改变A,B的持股比例，C股不变即可（本题没有用quadprog语句求解，最终结果与课本答案有偏差）。