1、数学实验实验九非线性规划大学数学实验作业非线性规划班级： 姓名： 学号： 日期： 目录【实验目的】 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_1形式命名，其中ex代表作业，9_4_1表示第九章第四小题第一个程序。
2、 自编函数均以exf9_4_1形式命名，exf代表作业函数，9_4_1表示第九章第四题第一个自编函数。【实验目的】1. 掌握用MATLAB优化工具箱和LINGO解非线性规划的方法；2. 练习建立实际问题的非线性规划模型。【实验内容】题目1（课本习题第九章第4题）某公司将3种不同含硫量的液体原料（分别记为甲、乙、丙）混合生产两种产品（分别记为A，B）。按照生产工艺的要求，原料甲、乙必须首先倒入混合池中混合，混合后的液体再分别与原料丙混合生产A，B。已知原料甲、乙、丙的含硫量分别是3%，1%，2%，进货价格分别为6千元/t，16千元/t，10千元/t；产品A，B的含硫量分别不能超过2.5%，1.5
3、%，售价分别为9千元/t，15千元/t。根据市场信息，原料甲、乙、丙的供应量都不能超过500t；产品A，B的最大市场需求量分别为100t，200t。(1) 应如何安排生产？(2) 如果产品A的最大市场需求量增长为600t，应如何安排生产？(3) 如果乙的进货价格下降为13千元/t，应如何安排生产？分别对（1）、（2）两种情况进行讨论。【第（1）问求解】【模型建立】(1)模型该题为带约束非线性规划问题，其模型包含决策变量、目标函数和约束条件。决策变量：设生产A所用的甲、乙、丙原料的量分别为；生产B所用的甲、乙、丙原料的量分别为；所生产的A、B产品的量分别为。单位均为t。目标函数：设总获利为z，同
4、时设销售收入为，进货费用为（单位均为千元）。则总获利由题中条件，销售收入进货费用故可得，总获利即为本题的目标函数。(2)约束条件浓度约束：按照生产工艺的要求，原料甲、乙必须先混合再分别与原料丙混合生产A，B。故有含硫量约束：要求产品A，B的含硫量分别不能超过2.5%，1.5%。故可得即原料供应约束：原料甲、乙、丙的供应量都不能超过500t。故得市场需求约束：产品A，B的最大市场需求量分别为100t，200t。故有A、B产量约束：生产的A、B产品的量与所用的甲、乙、丙原料的量所满足的关系非负约束：(3)模型整理将前面分析得到的带约束非线性规划模型化为如下标准形式：其中，将目标函数系数以及各约束条
5、件的数值代入标准形式，可得各函数表达式及参数数值分别为A1=1 0 0 1 0 0 0 0; 0 1 0 0 1 0 0 0; 0 0 1 0 0 1 0 0; 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 1; 3 1 2 0 0 0 -2.5 0; 0 0 0 3 1 2 0 -1.5b1=500 500 500 100 200 0 0A2=1 1 1 0 0 0 -1 0; 0 0 0 1 1 1 0 -1b2=0 0v1=0 0 0 0 0 0 0 0【模型求解及结果分析】根据已经建立好并化为标准形式的模型，利用Matlab优化工具箱对该模型进行求解。对于第（1）问，编
6、写运行代码如下：%-作业题9_4目标函数M文件（含梯度以及黑塞矩阵）exf9_4_1-function f,g,H=exf9_4_1(x)f=6*x(1)+16*x(2)+10*x(3)+6*x(4)+16*x(5)+10*x(6)-9*x(7)-15*x(8);% 给出分析梯度if nargout1 g=6,16,10,6,16,10,-9,-15;end% 给出黑塞矩阵（由于本题不能使用大规模算法，因此该黑塞矩阵不会被使用）if nargout2 H=zeros(8,8);end%-作业题9_4约束函数M文件(含梯度)exf9_4_2-function c1,c2,gc1,gc2=exf9
7、_4_2(x)c1=0;c2=x(1)*x(5)-x(2)*x(4);% 给出分析梯度if nargout2gc1=0 0 0 0 0 0 0 0;gc2=x(5),-x(4),0,-x(2),x(1),0,0,0 ;end%-作业题9_4函数M文件源程序ex9_4-clc;clear;% 输入不等式约束A1=1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 3 1 2 0 0 0 -2.5 0 0 0 0 3 1 2 0 -1.5;b1=500 500 500 100 200 0 0
8、;% 输入等式约束A2矩阵A2=1 1 1 0 0 0 -1 0 0 0 0 1 1 1 0 -1;b2=0 0; %等式约束% 迭代初值x0=10 10 10 10 10 10 10 10; % 自变量下限v1=0 0 0 0 0 0 0 0;% 采用中等规模算法模式opt=optimset(largescale,off);% 利用fmincon函数求解x,fv,ef,out,lag,grad,hess=fmincon(exf9_4_1,x0,A1,b1,A2,b2,v1,exf9_4_2,opt)得到运行结果为:最优解f=-400（最大利润z=-f=400，单位：千元）exitflag=1
9、 （收敛）out =iterations: 12（迭代次数）funcCount: 108（函数调用次数）algorithm: medium-scale: SQP, Quasi-Newton, line-search（所使用的算法为中等规模SQP方法）grad=6, 16, 10, 6, 16, 10-9, -15 （输出x点的梯度grad）hess=限于篇幅，数据略 （输出x点黑塞矩阵hess）从而可得，应进乙、丙各100t，且全部用于生产B产品200t，可获得最大利润40万元。【第（2）问求解】当产品A的最大市场需求量增长为600t时：则需要修改A产品市场需要需求为即b1向量变为b1=500
10、 500 500 600 200 0 0，并适当调整初值，比如调为x0=300 0 200 100 200 300 500 500。因此只需要在主程序中作相应修改，将不等式约束范围变化即可求解。得到结果如下：最优解f=-600 （最大利润z=-f=600，单位：千元）从而可得，当产品A的最大市场需求量增长为600t时，应该进货甲、丙各300t且全部用来生产A产品600t，可获得最大利润60万元.【第（3）问求解】如果乙的进货价格下降为13千元/t：则利润函数应变为f=9+156，只需修改目标函数M文件即可。若在第（1）问的条件下，即此时A的最大市场需求量为100t时，运行结果为最优解f=-75
11、0 （最大利润z=-f=750，单位：千元）从而可得，如果乙的进货价格下降为13千元/t且当产品A的最大市场需求量为100t时，应该进货甲50t、乙150t且全部用来生产B产品200t，可获得最大利润75万元。若在第（2）问的条件下，即此时A的最大市场需求量为600t时，运行结果为最优解f=-750 （最大利润z=-f=750，单位：千元）从而可得，如果乙的进货价格下降为13千元/t且当产品A的最大市场需求量为600t时，应该进货甲450t、乙150t且全部用来生产A产品600t，可获得最大利润75万元.综上所述，可得最优生产安排如下面得表1：表1：不同情形下的最优生产安排情况原料甲/t原料乙
12、/t原料丙/t产出A/t产出B/t总利润/万元第一问0100(全生产B)100(全生产B)020040第二问300(全生产A)0300(全生产A)600060第三问（1）50(全生产B)150(全生产B)0020075第三问（2）450(全生产A)150(全生产A)0600075【拓展实验、思考、对比、分析】1 代初值的对比、分析可以借助本题模型研究迭代初值对求解结果的影响。以本题第（1）问的模型为例，选取不同初值，得到结果如下面的表2：表2：第（1）问不同初值时的求解结果初值x0最优解x最大利润值/万元迭代次数（iterations）函数调用次数（funcCount）最优解处的梯度grad0
13、,0,000,0,0,0,0.6667,0.6667,0,1.33332.66673276, 16, 10, 6, 16, 10-9, -15100,1001000,0,0,0,100,100,0,20040121086, 16, 10, 6, 16, 10-9, -150,0,0,0,99,9,0,0,1990,0,0,0,100,100,0,200405456, 16, 10, 6, 16, 10-9, -15-100,-100-100接近零向量几乎为02186, 16, 10, 6, 16, 10-9, -15由表可以观察到如下结果：初值在一定范围内，得到的解是一致的。如表第3和第2行。
14、此时初值均在经历若干次迭代后最终收敛于同一个最优解。不同的初值迭代次数不一样，且解越靠近最优解迭代次数越少，如表第3行。这也是很容易理解的（选取初值不当错误停止的情况除外，比如第4行）。因此，在求解问题的过程中，应当尽量选取靠近最优解的点作为初值，从而提高程序运行效率。初值超过了一定范围后，可能无法得到符合实际要求的最优解甚至迭代过程不收敛。如表第1、4行。因此初值的选取会影响到所求得结果。可见对于某些算法而言，初值的选取是很重要的。梯度算法的对比、分析本题在求解过程中默认使用了目标函数和约束函数的数值梯度。也可以通过改变控制参量使用分析梯度。（也可以自定义目标函数的黑塞矩阵，但由于本题模型含
15、有非线性约束，故不能使用大规模算法，因而对于本题而言，在函数M文件中输入的黑塞矩阵不会被使用）。在目标函数M文件以及约束函数M文件中分别写入目标函数和约束函数的分析梯度（根据nargout的数值决定函数被调用时是否输出分析梯度），具体见前面所列的函数M文件代码。然后运行opt2=optimset(opt,GradObj,on,GradCon,on);x,fv,ef,out,lag,grad,hess=fmincon(exf8_10_1,x0,A1,b1,A2,b2,v1,exf8_10_2,opt2)得到的结果与数值梯度算法所得结果一起列表如下（以第一问中的模型为例进行对比分析）：表3：数值梯
16、度和分析梯度所得结果梯度算法最优解x最大利润值/万元迭代次数（iterations）函数调用次数（funcCount）数值梯度0,0,0,0,100,100,0,2004012108分析梯度0,0,0,1.027,101.027,97.946,0,20039.791427可见，数值梯度计算结果与分析梯度计算结果略有不同。迭代次数和函数调用次数也不同。但对于本题模型而言，数值梯度和分析梯度差别不大。对拉格朗日乘子的探索、分析可以借助本题模型，仿照前一题对结果输出的lag结构变量即程序停止时的相应约束的拉格朗日乘子作一些分析。拉格朗日乘子包含了程序运行的有关信息，其维数等于约束条件的个数，其非零分
17、量对应于起作用的约束，即等号严格成立时的约束。比如，以第一问模型为例，运行lag.ineqlin可得到lag.ineqlin=0, 0, 0, 0, 2, 0, 6该乘子对应于约束不等式约束。此时如果将b1的第6个分量加1，即令b1=500 500 500 100 200 0 1，运行程序，可得到目标函数的值变为406。而改动前目标函数的值为400。406-400=6，该值恰为lag.Ineqlin第6分量的值。可见，拉格朗日乘子表示的是对应约束的右端项增加一个单位时（其它条件不变），目标函数的增加量（在一定的适用范围以内）。【本题小结】本题对应的模型为带约束非线性规划问题。利用Matlab的
18、优化工具箱可以求解本题的模型。求得：对于第1问：应进乙、丙各100t，且全部用于生产B产品200t，可获得最大利润40万元。对于第2问：应该进货甲、丙各300t且全部用来生产A产品600t，可获得最大利润60万元；对于第3（1）问：应该进货甲50t、乙150t且全部用来生产B产品200t，可获得最大利润75万元；对于第3（2）问：应该进货甲450t、乙150t且全部用来生产A产品600t，可获得最大利润75万元。可以借助本题模型分析初值对求解结果的影响。在一定范围内，不同初值得到的解是一致的。但不同的初值迭代次数不一样，且解越靠近最优解迭代次数越少；初值超过了一定范围后，可能无法得到符合实际要
19、求的最优解甚至迭代过程不收敛。因此，在求解问题的过程中，应当尽量选取靠近最优解的点作为初值，从而提高程序运行效率。本题可以使用目标函数以及约束函数的数值梯度和分析梯度分别求解。数值梯度计算结果与分析梯度计算结果完全一致，但迭代次数和函数调用次数不同。可以借助本题模型继续对拉格朗日算子的含义作一些分析。其零分量表示与之对应的约束条件不起作用，而非零分量对应于等号严格成立时的约束，且有一定的实际含义：拉格朗日乘子表示的是对应约束的右端项增加一个单位时（其它条件不变），目标函数的增加量（当然，仅在一定的适用范围以内成立）。题目2（课本习题第九章第8题）美国某3种股票（A, B, C）12年（1943
20、1954年）的价格（已经包括了分红在内）每年的增长情况如表1所示（表中还给出了相应年份的500种股票的价格指数的增长情况）。例如，表中第一个数据1.300的含义是股票A在1943年的年末价值是其年初价值的1.300倍，即收益为30%，其余数据的含义依此类推。假设你在1955年时有一笔资金准备投资这三种股票，并期望年收益率至少达到15%，那么你应当如何投资？此外，考虑以下问题：（1）当期望的年收益率在10%100%之间变化时，投资组合和相应的风险如何变化？（2）假设除了上述三种股票外，投资人还有一种无风险的投资方式，如购买国库券。假设国库券的年收益率为5%，如何考虑该投资问题？（3）假设你手上目
21、前握有的股票比例为：股票A占50%，B占35%，C占15%。这个比例与你得到的最优解可能有所不同，但实际股票市场上每次股票买卖通常总有交易费，例如按交易额的1%收取交易费，这时你是否仍需要对手上的股票进行买卖（换手），以便满足“最优解”的要求？表4年份股票A股票B股票C股票指数19431.3001.2251.1491.25899719441.1031.2901.2601.19752619451.2161.2161.4191.36436119460.9540.7280.9220.91928719470.9291.1441.1691.05708019481.0561.1070.9651.05501
22、219491.0381.3211.1331.18792519501.0891.3051.7321.31713019511.0901.1951.0211.24016419521.0831.3901.1311.18367519531.0350.9281.0060.99010819541.1761.7151.9081.526236【模型建立】 设购买A, B, C三种股票的比例分别为，记股票A, B, C的收益率分别为，根据题意，是一个随机变量，投资的总收益率也是一个随机变量。用和分别表示随机变量的数学期望和方差，表示两个随机变量的协方差，则根据1943年1954年的统计数据，用excel求解期望、
23、方差、协方差可以得到以下结果：（注：使用Excel软件计算方差有两种函数：VAR函数和VARP函数，两种方法所用公式不同，适用场合也不同：VAR函数的计算公式为：，代表样本方差；VARP函数的计算公式为：，代表总体方差。而在本实验条件下，我选用的是VARP函数）表5：excel求解各股票期望、方差及协方差年份股票A股票B股票C19431.30.31.2250.2251.1490.14919441.1030.1031.290.291.260.2619451.2160.2161.2160.2161.4190.41919460.954-0.0460.728-0.2720.922-0.07819470
24、.929-0.0711.1440.1441.1690.16919481.0560.0561.1070.1070.965-0.03519491.0380.0381.3210.3211.1330.13319501.0890.0891.3050.3051.7320.73219511.090.091.1950.1951.0210.02119521.0830.0831.390.391.1310.13119531.0350.0350.928-0.0721.0060.00619541.1760.1761.7150.7151.9080.908E1.089083330.089083331.213666670.2
25、13666671.234583330.23458333D0.009906910.009906910.053525720.053525720.086374580.08637458covcovAB 0.011373covBC 0.050807covAC 0.011985于是，投资的总期望收益率为.用总收益率的方差衡量投资的风险。投资总收益率的方差为根据题目有以下约束条件：购买比例约束： 非负约束： 一般来说，希望投资的收益率高而风险小，即最大而最小，是一个多目标优化问题，通常需转化为单目标优化问题来求解。要求年收益率至少为15%，则有：因此模型可以化为： s.t. 为了利用MATLAB优化工具箱求
26、解二次规划，所以将模型化为如下形式：s.t. 其中，。【模型求解】根据上述模型编写程序如下：%-作业题9_8函数M文件源程序ex9_8_1-clear;clc;H=0.0198 0.0227 0.0240 0.0227 0.1070 0.1016 0.0240 0.1016 0.1728; %给出目标函数二次型的矩阵A1=-0.0891 -0.2137 -0.2346;b1=-0.15;A2=ones(1,3);b2=1;v1=zeros(1,3); %约束条件x0=0.5,0.3,0.1; %给出初值x,fv,ef,out,lag=quadprog(H,0 0 0,A1,b1,A2,b2,v
27、1,x0)输出结果：运行程序输出结果如下：x = 0.5301 0.3573 0.1126fv = 0.0205 即A, B, C三种股票的购买比例分别为53.01%、35.73%、11.26%时，投资风险最小，方差为0.0205。【第（1）问求解】求当期望的年收益率在10%100%之间变化时投资组合和相应的风险的变化情况。改变约束条件，用for循环来实现年收益率的修改，b1从10%开始以0.1%的间隔增加到100%，输出每次改变后风险变化的情况及相应的投资组合情况。修改模型如下s.t. Matlab程序如下：%-作业题9_8函数M文件源程序ex9_8_2-clear;clc;H=0.0198
28、 0.0227 0.0240 0.0227 0.1070 0.1016 0.0240 0.1016 0.1728; A1=-0.0891 -0.2137 -0.2346;A2=ones(1,3);b2=1;v2=ones(1,3);v1=zeros(1,3);x0=0 0 0;for i =1:1:901b1=-0.001*i-0.099;x(:,i),fv(i)=quadprog(H,0 0 0,A1,b1,A2,b2,v1,v2,x0)endy=0.1:0.001:1;figure(1)plot(y,x(1,:),-),gridxlabel(期望收益率)ylabel(最佳A股比例)title(期望收益率10%100%变化时，最小投资
