书中Matlab源程序 第1 章 绪论
【例1-1】有一名学生，期末有5门功课要考试，可用的复习时间有18小时。假定这五门课程分别是数学、英语、计算机基础、画法几何和专业概论。如果不复习直接参加考试，这五门功课预期的考试成绩分别为65分、60分、70分、60分和65分。复习以1小时为一单元，每增加1小时复习时间，各门功课考试成绩就有可能提高，每复习1小时各门功课考试成绩提高的分数分别为3分、4分、5分、4分和6分。问如何安排各门功课的复习时间可使平均成绩不低于80分，并且数学和英语成绩分别不低于70分和75分。
解：设分配在数学、英语、计算机基础、画法几何和专业概论这五门功课的复习时间分别为x1,x2,x3,x4,x5，则可列出如下的目标函数和限制条件为：
f(x)?x1?x2?x3?x4?x5 x1?x2?x3?x4?x5?18
(3x1?4x2?5x3?4x4?6x5?320)/5?80
3x1?65?70
本例具体程序如下：
%li_1_1 f=[1 1 1 1 1];
A=[1 1 1 1 1; -3 -4 -5 -4 -6; -3 0 0 0 0; 0 -4 0 0 0; 3 0 0 0 0; 0 4 0 0 0; 0 0 5 0 0; 0 0 0 4 0; 0 0 0 0 6]; b=[18;-80;-5;-15;35;40;30;40;35]; lb=zeros(4,1)
[x,fval]=linprog(f,A,b,[],[],lb)
计算结果为：
x = 1.6667 3.7500 5.0000 0.0000 5.8333 fval =
16.2500
1
【例1-2】某工厂要生产两种规格的电冰箱，分别用Ⅰ和Ⅱ表示。生产电冰箱需要两种原材料A和B，另外需设备C。生产两种电冰箱所需原材料、设备台时、资源供给量及两种产品可获得的利润如表1-1所示。问工厂应分别生产Ⅰ、Ⅱ型电冰箱多台，才能使工厂获利最多？
表1.1 资源需求与限制
资源 设备 原料A 原料B 单位产品获利 产品Ⅰ用原料限制 Ⅰ 1 2 0 220元 Ⅱ 1 1 1 250元 资源限制 1200台时 1800千克 1000千克 求最大收益 ?800千克 解：设生产Ⅰ、Ⅱ两种产品的数量分别为x1,x2。则可获得的最大收益为
maxf(x)? 220 x1 ? 250 x2 ,x?R2
s.. tx1 ? x2 ? 1200 2 x1 ? x2 ? 1800 x1 ? 800x2 ? 1000 x1 , x2 ? 0 Matlab求解程序如下：
%li_1_2 clc; close all; f=-[220 250]; A=[1 1;2 1;1 0;0 1]; b=[1200;1800;800;1000]; xl=[0 0];
[x,fval]=linprog(f,A,b,[ ],[ ],xl) x1=[0:1800]; x2=[0:2000];
[xm1,xm2]=meshgrid(x1,x2); x21=1200-x1; x22=1800-2*x1; x23=(-fval-220*x1)/250;
2
plot(x1,x21,x1,x22,[0:1:1000],1000,800,[0:1:1500],x1,x23,'r') axis([0,1400,0,2000]) xlabel('x1'); ylabel('x2'); hold on
z=200*xm1+250*xm2; [C,h]=contour(xm1,xm2,z); text_handle = clabel(C,h);
set(text_handle,'BackgroundColor',[1 1 .6],'Edgecolor',[.7 .7 .7]); hold off
【例1-3】绘制下面函数的曲线。
y(x)?2sin(x)?ln(x)
解：应用plot()函数绘制该函数曲线的程序如下：
%li_1_3
f=inline('2*sin(x)+log(x)','x') x=linspace(0.1,2*pi,15); y=feval(f,x);
plot(x,y,'-rs','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g', 'MarkerSize',10) xlabel('0.1\\leq \\Theta \\leq 2\\pi') ylabel('2sin(\\Theta)+ln(\\Theta)'); title('Plot of 2sin(\\Theta)+ln(\\Theta)')
text(pi/4,sin(-pi/4),'\\leftarrow 2sin(\\Theta)+ln(\\Theta)','HorizontalAlignment','left') legend('-') grid on
【例1-4】用图形表示如下优化模型，并求解。 minf(x)?4x1?x2?12
2s.t.x12?xx?2532(x1?5)?(x2?5)?16?022
解：该绘制目标函数曲面、约束函数曲线及求解程序如下：
（1）绘制目标函数曲面的程序
%li_1_4_1 function li_1_4_1() clc; clear all; close all; n=20;
3
x1=linspace(0,2,n); x2=linspace(0,6,n);
[xm1,xm2]=meshgrid(x1,x2); for i=1:n for j=1:n
xx=[xm1(i,j),xm2(i,j)]; zm(i,j)=fun_obj(xx); end end figure(1)
meshc(xm1,xm2,zm) xlabel('x1'); ylabel('x2'); zlabel('zm')
（2）绘制约束函数曲线及求解的程序
%li_1_4_2 function li_1_4_2() clc; clear all; close all; x0=[1,1];
[x,fval,exitflag,output]=fmincon(@fun_obj,x0,[ ],[ ],[ ],[ ],[ ],[ ],@fun_cons) n=20;
x1=linspace(0,6,n); x2=linspace(0,10,n); [xm1,xm2]=meshgrid(x1,x2); for i=1:n for j=1:n
xx=[xm1(i,j),xm2(i,j)]; zm(i,j)=fun_obj(xx); end end figure(1)
f1=inline('sqrt(25-x.^2)','x'); f2=inline('sqrt(16-(x-5).^2)+5','x'); y1=feval(f1,x1); y2=feval(f2,x1);
y3=sqrt((4*x1.^3)-12-fval+0.01)
4
plot(x1,y1,x1,y2); hold on
plot(x1(1:8),y3(1:8),'--r') hold on
[c,h]=contour(xm1,xm2,zm,20); clabel(c,h); xlabel('x1'); ylabel('x2');
function f=fun_obj(x) f=4*x(1)^3-x(2)^2-12; function [c,ceq]=fun_cons(x)
c=[x(1)^2+x(2)^2-10*x(1)-10*x(2)+34 -x(1) -x(2)];
ceq=[x(1)^2+x(2)^2-25];
5
第2章 优化设计的数学基础
【例2-11】 用图解法分析下面面优化问题的凸性，并求其最小值。
2minf(x)?(x1?3)2?x2s..tg1(x)?x12?x2?4?0 （2-28）
g2(x)??x1?0g3(x)??x2?0解：根据所给目标函数和约束函数函数，编制如下程序：
function kt_test1_plot1 clc; clear all; close all;
x=linspace(0,2.5,30); [xm,zm1]=meshgrid(x,[0,6]); y=4-x.^2;
ym=meshgrid(y,[0,20]); mesh(xm,ym,zm1); hold on r=2;
t=linspace(0,2*pi,45); x=r*cos(t)+3; y=r*sin(t);
[xm,ym]=meshgrid(x,y); zm2=(xm-3).^2+ym.^2; mesh(xm,ym,zm2); hold on
xx=linspace(-1,6,30); yy=linspace(-2,5,30);
[xxm,zzm]=meshgrid(0*xx,[0,6]); [yym]=meshgrid(yy,[-10,0]); mesh(xxm,yym,zzm); hold on
[yym,zzm]=meshgrid(0*yy,[0,6]); [xxm]=meshgrid(xx,[-10,0]); mesh(xxm,yym,zzm); axis([-1,6,-3,5,-2,12])
6
xlabel('x');ylabel('y');zlabel('z'); hold off figure(2)
x=linspace(0,6,30); y=linspace(0,5,30); [xm,ym]=meshgrid(x,y); zm3=(xm-3).^2+ym.^2; y1=4-x.^2; plot(x,y1,'k'); hold on
[c,h]=contour(xm,ym,zm3,10); text=clabel(c); hold off
xlabel('x');ylabel('y'); axis([0 6 0 5 ]);
【例2-12】分析式（2-29）和式（2-30）所示非线性有约束最小值问题。
minf(x)??(x1?1)2?(x2?10)2?102s..tg1(x)?x12?x2?1?0?x1?0?x2?0 （2-29）
minf(x)?(x1?1)2?(x2?10)2?102s..tg1(x)?x12?x2?1?0?x1?0?x2?0 （2-30）
解：首先绘制目标函数和约束函数曲面和曲线。以式（2-29）为例，绘制函数曲面和曲线的程序如下：
function kt_test1_plot2 clc;
thita=linspace(0,2*pi,50); r=1;
x=r*cos(thita); y=r*sin(thita);
[xm,zm1]=meshgrid(x,[-5,10]); [ym]=meshgrid(y,[0,20]); figure(1);
7
mesh(xm,ym,zm1); hold on
[xm,ym]=meshgrid(linspace(-2,2,30),linspace(-2,2,30)); zm2=(-(xm-1).^2-(ym-2).^2+10); mesh(xm,ym,zm2); hold on
axis([-2,2,-2,2,-15,12]) xlabel('x');ylabel('y');zlabel('z'); hold off figure(2)
x=linspace(-1,1,30); y1=sqrt(1-x.^2); y2=-sqrt(1-x.^2); plot(x,y1,'k',x,y2,'-r'); hold on
[c,h]=contour(xm,ym,zm2,20); text=clabel(c); hold on yy=[-2:0.1:2]
plot(zeros(length(yy)),yy,'--'); hold on
plot(yy,zeros(length(yy)),'--') hold off
xlabel('x');ylabel('y'); axis([-2 2 -2 2 ]);
8
第 3 章 线性规划
【例3-1】用单纯形法求解下面的线性规划问题。
maxf(x)?2x1?x2?2x3s.t.x1?x2?2x3?5x1?3x2?x3?3x1,x2,x3?0解：
（1）先用MATLAB线性规划函数求解，其计算程序如下：
%ch31_li1 clc; close all; f=-[2 1 -2]; A=[1 1 2;1 3 -1]; b=[5 3]; xl=[0 0 0];
[x,fval]=linprog(f,A,b,[ ],[ ],xl)
3.3 单纯形法的Matlab程序及实例
程序清单如下：
function [x,fmax]=linear_pro_max(cf,cb,A,b,indexb1) n=length(cf); max_sigma=1; m=length(cb); indexb=indexb1; theta=zeros(size(m,1)); while max_sigma>0 for j=1:n
sigma(j)=cf(j)-sum(cb(:).*A(:,j));
end
max_sigma=max(sigma); if(max_sigma>0)
pvj=find(sigma==max_sigma); theta=b./A(:,pvj); min_theta=min(theta); max_sigma min_theta
9
pvi=find(theta==min_theta); cb(pvi)=cf(pvj); indexb(pvi)=pvj; pvi pvj cb cf for i=1:m if(i~=pvi) for j=1:n
AA(i,j)=A(i,j)-A(i,pvj)*A(pvi,j)/A(pvi,pvj); end
bb(i)=b(i)-A(i,pvj)*b(pvi)/A(pvi,pvj); else
AA(i,:)=A(i,:)/A(pvi,pvj); bb(i)=b(i)/A(pvi,pvj); end end end
A=AA;b=bb'; end s=1:n; x=zeros(n,1); for i=1:m
k=find(s==indexb(i)); if(k~=0) x(k)=b(i); end end fmax=cf*x;
【例3-4】用单纯形法Matlab程序求解例1-2。
解：为方便起见，重新列出所给问题线性规划的标准形：
maxf(x)? 220 x1 ? 250 x2 ?0x3?0x4?0x5?0x6
s.. tx1 ? x2 +x3? 1200 2 x1 ? x2 +x4? 1800 x1 ?x5= 800x2 +x6= 1000
10
xj?0,j?1,2,?6
编写用户程序。cf为目标函数中变量系数行向量；cb为初选基变量行向量；indexb1为初始基变量下标索引；A为等式约束方程系数矩阵；b为等式约束方程右端项。应户程序为：
function linear_pro_max_test1 clc clear all;
cf=[220 250 0 0 0 0]; cb=[cf(3),cf(4),cf(5),cf(6) ]; indexb1=[3 4 5 6]; A=[1 1 1 0 0 0 2 1 0 1 0 0 1 0 0 0 1 0
0 1 0 0 0 1]; b=[1200;1800;800;1000];
[x,fmax]=linear_pro_max(cf,cb,A,b,indexb1)
计算结果：
x=[200 1000]。
3.5改进的单纯形法的Matlab程序及实例
根据修正的单纯形法计算步骤，应用Matlab语言编写计算程序，程序清单如下。
function [x fmax]=linear_promax_rev(cf,cb,A,b,indexa1,indexb1) cf1=cf; n=length(cf); m=length(cb); indexa=indexa1; indexb=indexb1; x0=b; for j=1:n A1(j)={A(:,j)}; end
e0=[A1{m:n}]; e0inv=inv(e0); xe0=e0inv*x0; r0=1; kk=0; while r0>0
11
fprintf('===========================================================') kk=kk+1 indexa,indexb cb e0inv for i=1:n-m k=indexa(i); cc=A1{k};
r(i)=cf(k)-cb*e0inv*cc; end r0=min(r); max_sigma=max(r); pvj=find(r==max_sigma); rr=e0inv*A1{pvj}; theta=x0./rr;
min_theta=min(theta); pvi=find(theta==min_theta); temp=cb(pvi); cb(pvi)=cf(pvj); cf(pvj)=temp;
indexa(pvj)=indexb(pvi); indexb(pvi)=pvj; e=A1{pvj}; e1=e0; for i=1:m e1(i,pvi)=e(i); end
e1inv=inv(e1); xe1=e1inv*b; e0=e1; e0inv=e1inv; x0=xe1;
max_sigma,min_theta,pvi,pvj e1 e1inv
fprintf('===========================================================') end s=1:n; x=zeros(n,1);
12
for i=1:m
k=find(s==indexb(i)); if(k~=0) x(k)=xe1(i); end end fmax=cf1*x;
【例3-6】 应用单纯形法Matlab程序计算例3-5。
解：将线性规划问题写成标准型，编写用户程序。cf为目标函数中变量系数行向量；
cb为初选基变量行向量；indexa1为初始非基变量下标索引；indexb1为初始基变量下标索
引；A为等式约束方程系数矩阵；b等式约束方程右端项。
用户程序如下：
function linear_promax_rev_test2 clc; clear all; cf=[50 100 0 0 0]; A=[1 1 1 0 0 2 1 0 1 0 0 1 0 0 1];
b=[300;400;250]; cb=[cf(3),cf(4) cf(5)]; indexa1=[1 2]; indexb1=[3 4 5];
[x fmax]=linear_promax_rev(cf,cb,A,b,indexa1,indexb1)
计算结果为：
x=[50 250 0 50 0],fmax=27500。
13
第 4 章 一维搜索方法
2. 进退法的MATLAB程序
该程序是按照寻找最佳步长编写的，其数学模型为式（4-2）。 程序说明如下：
函数opt_range_serach(xk0,dir0,th)中输入参数； xk0：初始点；
dir0：给定的搜索方向； th：步长增量。
%opt_range_serach1.m
function [opt_step,fo,xx] = opt_range_serach1(f,xk0,dir0,th)
%用进退法搜索三个点，使中点函数值最小；输出步长，函数值，设计变量值 %xk0:初始点 %th:步长 t1=0;t2=th;
xk1=xk0;xk2=xk1+t2*dir0; x0=xk1; f1=feval(f,x0); x0=xk2; f2=feval(f,x0); if f2 th=-2*th;t3=t1;f3=f1;t1=t2;f1=f2;t2=t3;f2=f3; t3=th; xk3=xk1+t3*dir0; x0=xk3; f3=feval(f,x0); end ii=0; while f2>f3 t1=t2;f1=f2;t2=t3;f2=f3; t3=t2+th; xk3=xk1+t3*dir0; x0=xk3; 14 f3=feval(f,x0); end xx1=xk1+t1*dir0; xx2=xk1+t2*dir0; xx3=xk1+t3*dir0; if th<0 opt_step=[t3 t2 t1]; xx=[xx3 xx2 xx1]; fo=[f3 f2 f1]; else opt_step=[t1 t2 t3]; xx=[xx1 xx2 xx3]; fo=[f1 f2 f3]; end end 【例4-1】用进退法计算函数 f(x)?2?x2 的单峰区间，初始点x0?2。 解：用户程序如下： %range_search_test1 clc; f=inline('2+x^2','x'); xk0=2;th=0.5;dir0=1; [opt_step,fo,xx] = opt_range_serach1(f,xk0,dir0,th) opt_step = 结果为： opt_step=[-3 -2 -1] xo=[-1 0 1] fo=[3 2 3] 1. 黄金分割法的计算框图 图4.5是黄金分割法的计算框图。图中?1,?2为给定的任意小的精度，k为区间缩短的次数。 2. MATLAB程序 %golden_search function [xo,fo] =golden_search(f,a,b,r,TolX,TolFun,k) kk=1; 15 while kk>0 h = b - a; rh = r*h; c = b - rh; d = a + rh; fc = feval(f,c); fd = feval(f,d); if k <= 0 || abs(h) < TolX && abs(fc - fd) xo = d;fo = fd; kk=0; end if k == 0;fprintf('达到计算次数');kk=0; end else if fc < fd b=d;k=k-1; else a=c;k=k-1; end end end 【例4-2】计算目标函数 f(x)?2?x2 在区间[?2,2]内的极小点。 解：用户程序如下： %golden_s_test1.m function golden_s_test1 clc; clear all; gs_fun = inline('2+x^2','x'); a =-2; b = 2; r =(sqrt(5)-1)/2; TolX =1e-7; TolFun = 1e-7; MaxIter =50; [xo,fo]=fminbnd(gs_fun,a,b) [xo,fo] =golden_search(gs_fun,a,b,r,TolX,TolFun,MaxIter) 计算结果： xo= 1.4142e-0082; fo=2 2. 二次函数插值法的迭代过程与程序框图 16 MATLAB程序如下： function [opt_step,fo,xo]=opt_step_quad2(f,xk0,dir0,th,TolX,TolFun,MaxIter) %opt_step_quad.m [t012,fo,xx] = opt_range_serach1(f,xk0,dir0,th); %search for the optimum step corresponding to minimum f(x) by quadratic approximation method k=MaxIter; while k>0 k=k-1; t0=t012(1);t1=t012(2);t2=t012(3); x0 = xk0+t012(1)*dir0; x1 =xk0+ t012(2)*dir0; x2 =xk0+t012(3)*dir0; f0 =feval(f,x0); f1 =feval(f,x1); f2 =feval(f,x2); nd = [f0 - f2 f1 - f0 f2 - f1]*[t1*t1 t2*t2 t0*t0; t1 t2 t0]'; t3 = nd(1)/2/nd(2); x3=xk0+t3*dir0;f3 = feval(f,x3); if k <= 0 | abs(t3 - t1) < TolX | abs(f3 - f1) < TolFun opt_step=t3; xo=xk0+opt_step*dir0; fo=f3; return else if t3 < t1 if f3 < f1, t012 = [t0 t3 t1]; f012 = [f0 f3 f1]; else t012 = [t3 t1 t2]; f012 = [f3 f1 f2]; end else if f3 <= f1, t012 = [t1 t3 t2]; f012 = [f1 f3 f2]; else t012 = [t0 t1 t3]; f012 = [f0 f1 f3]; end end end end opt_step=t3; xo=xk0+opt_step*dir0; fo=f3; end 【例4-3】用二次插值法计算目标函数 f(x)?2?x2 在区间[?2,2]内的极小点及最佳步长。 17 解：用户程序如下： %opt_step_quad_test1 clc; f=inline('2+x^2','x'); xk0=2;th=0.5;dir0=1;TolX=1e-4;TolFun=1e-4;MaxIter=50; [opt_step,fo,xx]=opt_step_quad2(f,xk0,dir0,th,TolX,TolFun,MaxIter) 结果为： opt_step=[-2] xo=[0] fo=[2] 【例4-5】根据Lagrange插值公式，用MATLAB编程计算ln(x)在x?1.5的的值，差值区间x?[1,20]。 %int_lagrange function int_lagrange clc; x=1:0.5:20; y=log(x); xx=1.5; n=length(x); ff=0; for i=1:n p=1; for j=1:n if j~=i p=p*(xx-x(j))/(x(i)-x(j)); end end ff=ff+p*y(i); end [xx,ff] log(1.5) 计算结果为 x=1.5000,f(x)=0.4055。 3）牛顿迭代法的MATLAB程序 function [x,fx] = root_newton(f,x0,TolX,MaxIter) %root_newton.m 18 TolFun=eps; xx = x0; fx = feval(f,x0) for k = 1: MaxIter dfdx =fun_dff1(f,xx) %numerical drv dx = -fx/dfdx; xx = xx+dx; fx = feval(f,xx); if norm(fx) if k == MaxIter, fprintf('达到最大迭代次数’), end function g=fun_dff1(f,x,dh) if nargin < 3, dh=1e-4;end h2 = 2*dh; M=length(f);N = length(x(1,:));I = eye(N); for n = 1:N g(:,n) = (feval(f,x + I(n,:)*dh)-feval(f,x - I(n,:)*dh))/h2; end 改进的牛顿迭代法的MATLAB程序如下。 function [xk1,fk1] =root_broyden(f,x0,TolX,MaxIter) %root_broyden.m to minimize an objective ftn f(x) by the broyden method. N=length(x0(1,:)); TolFun=TolX; ak=eye(N); n=0; while n<=MaxIter n=n+1; fk=feval(f,x0); xk1=x0-fk/ak; sk=xk1-x0; fk1=feval(f,xk1); yk=fk1-fk; deltak=(yk'-ak*sk')*sk/(sk*sk'+1e-7); ak1=ak+deltak; ak=ak1; x0=xk1; if norm(fk1) if n == MaxIter, fprintf('The best in %d iterations\\n',MaxIter), end 19 end 2）对分法的MATLAB程序 对分法的MATLAB程序如下。 function [alpha,fc]=root_bisect(f,a,b,x0,dir,TolFun) % root_bisect.m a1=x0+a*dir;b1=x0+b*dir; fa=feval(f,a1);fb=feval(f,b1); while norm(fb-fa)>TolFun fa=feval(f,a1);fb=feval(f,b1); c1=(a1+b1)/2; fc=feval(f,c1); if fa*fc<0 b1=c1; else a1=c1; end end alpha=(c1(1)-x0(1))/dir(1); 20 c=c+(g(i)<0)*log(abs(g(i))); end fc=f+r*c; 计算结果为： xo = fc = 1.00014459741663 1.00010436403355 1.99984323680675 例8-3的Matlab程序如下。 function pen_ex_test %外点罚函数法 clc; global r format long; x0=[0.5 0.5]; r=1e-6; c=1; k=0 TolX=1e-3 while c>1e-6 k=k+1; r=r*1.1 [x fc]=fminsearch(@fpen_in,x0) [fc f c]=fpen_in_sjm1(x) x0=x; end function [fc f c]=fpen_in_sjm1(x) global r f=x(1)^2+2*x(2)^2; c=-x(1)-x(2)+1; fc=f+r*(c>0)*c^2; 计算结果如下： x = fc = 0.66666599488580 0.66679835346354 0.33320078985086 【例8-4】用混合罚函数法求下面的优化问题。 46 2minf(x)?4x1?x2?12 2 s.t. g1(x)?-10x1?x1-10x2?x22?20?0 g2(x)??x1?0 g3(x)??x2?0 22 h1(x)?x1?x2?25?0 解：构造罚函数为 2?(x)?4x1?x2?122?r(1/(-10x1?x1-10x2?x22?20)?1/x1?1/x2) 2?1/r(x12?x2?25)2MATLAB程序如下： function pen_hy_fmtest2 clc; clear all; close all; format long; x0=[2,3]; r=1e8; df=1e10;; k=0; while df>0.00001 k=k+1; r=r*0.999; [x fc]=fminsearch(@fpen_hyfun,x0,[ ],r); [fc f c]=fpen_hyfun(x,r); df=abs(f-fc) x0=x; end k x,f function [fc f c]=fpen_hyfun(x,r) f=4*x(1)-x(2)^2-12; 47 c=[-10*x(1)+x(1)^2-10*x(2)+x(2)^2+20 -x(1) -x(2)]; ceq=25-x(1)^2-x(2)^2; fc=0; for i=1:3 fc=fc+(c(i)<0)*(1/c(i)); end fc=f-r*fc+1/sqrt(r)*(ceq)^2; 计算结果为：x=[4.158?10 5.000]，f = -36.998 -4 48 第9 章 多目标函数优化设计 【例9-1】 现有现金70元，可用来可用来购买菠萝和苹果。菠萝5元/kg，苹果3元/kg，要求总斤数不少于15kg，菠萝不少于5kg。 问：(1) 购买菠萝和苹果各多少斤，才能在满足要求的条件下花钱最少？(2) 购买菠萝和苹果各多少斤，才能在满足要求的条件下所买的菠萝和苹果最多？ 解：通俗地说，这是一个如何安排资金，少花钱多办事的问题。设购买菠萝x1kg，苹果x2kg。可以列出如下的优化数学模型： minf1(x)?5x1?3x2，x?R maxf2(x)?x1?x2，x?R 22s..t5x1?3x2?70 x1?x2?15 x1?5 x2?0 (1)如果只考虑目标函数minf1(x)?5x1?3x2，则应用Matlab求解的程序为： %li9_1 f=[5 3]; A=[5 3;-1 -1;-1 0]; b=[70;-15;-5]; xl=[0;0]; [x,f,exitflag]=linprog(f,A,b,[ ],[ ],xl) 最优解为：x=[5 10], f1(x)?55元，f2(x)?15kg。 (2)如果只考虑目标函数maxf2(x)?x1?x2，则应用Matlab求解的程序为： f=-[1 1]; A=[5 3;-1 -1;-1 0]; b=[70;-15;-5]; xl=[0;0]; [x,f,exitflag]=linprog(f,A,b,[ ],[ ],xl) 最优解为：x=[5 15], f1(x)?70元，f2(x)?20kg。 【例9-5】用线性加权法求解例9-1。 49 解：取两组权系数：1）w1?0.5,w2?0.5；2）w1?0.2,w2?0.8。相应的MATLAB计算程序如下： %li9_5 function li9_5 clc; clear all; A=[5 3;-1 -1;-1 0]; b=[70;-15;-5]; xl=[0,0];x0=[0,0]; w(1,1)=0.5;w(1,2)=0.5; w(2,1)=0.2;w(2,2)=0.8; for i=1:2 i [x,f,exitflag]=fmincon(@(x)fun_obj(x,w(i,:)),x0,A,b,[ ],[ ],xl) f1=5*x(1)+3*x(2) f2=x(1)+x(2) end f1=inline('(70-5*x)/3','x'); f2=inline('15-x','x'); x=0:0.1:20; plot(x,f1(x),x,f2(x),5*ones(length(x)),0:0.1:20) axis([0,16,0,25]) function f1f2=fun_obj(x,ww) f1f2=ww(1)*(5*x(1)+3*x(2))-ww(2)*(x(1)+x(2)); 50