1、MATLAB 简介数学建模 第六部分 MATLAB 优化算法一、线性规划算法调用格式：x, fval, exitflag= linprog(f,A,b, Aeq,beq,lb,ub, x0)说明：返回值 x为最优解向量，fval 为最优值；若没有不等式约束，则令 A= 、b= ；lb ,ub为变量 x的下界和上界，x0 为初值点； exitflag 描述函数计算的退出条件：若为正值，表示目标函数收敛于解 x处；若为负值，表示目标函数不收敛；若为零值，表示已经达到函数评价或迭代的最大次数。例 1、求解线性规划问题max f=70x1+120x2s.t 9x1+4x236004x1+5x22000
2、3x1+10x23000x1,x20将其转换为标准形式：min f=-70x1-120x2s.t 9x1+4x236004x1+5x220003x1+10x23000x1,x20算法如下：f=-70 -120;A=9 4 ;4 5;3 10 ;b=3600;2000;3000;lb=0 0;ub=;x,fval,exitflag=linprog(f,A,b,lb,ub)MATLAB 简介数学建模 maxf=-fval例 2、求解线性规划问题max f=0.15x1+0.1x2+0.08 x3+0.12 x4s.t x1-x2- x3- x40x2+ x3- x40x1+x2+x3+ x4=1x
3、j0 ， j=1,2,3,4将其转换为标准形式：min z=-0.15x1-0.1x2-0.08 x3-0.12 x4s.t x1-x2- x3- x40-x2- x3+ x40x1+x2+x3+ x4=1xj0 ， j=1,2,3,4算法如下：f = -0.15;-0.1;-0.08;-0.12;A = 1 -1 -1 -1；0 -1 -1 1; b = 0; 0;Aeq=1 1 1 1; beq=1; lb = zeros(4,1);x,fval,exitflag = linprog(f,A,b,Aeq,beq,lb)f=-fval二、二次规划算法调用格式： x,fval,exitflag
4、= quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)说明：输入参数中，x0 为初始点；若无等式约束或无不等式约束，就将相应的矩阵和向量设置为空。输出参数中，x 是返回最优解；fval 是返回解所对应的目标函数值；exitflag 是描述搜索是否收敛。 例 3、求解二次规划问题min f(x)= x1-3x2+3x12+4x22-2x1x2s.t 2x1+x22MATLAB 简介数学建模 -x1+4x23算法如下： f=1;-3;H=6 -2;-2 8;A=2 1;-1 4;b=2;3;X,fval,exitflag=quadprog(H,f,A,b)例 4、求解二次规划问题
5、min x12+2x22-2x1x2-4x1-12x2s.t x1+x22-x1+2x22 2x1+x23x10, x 20算法如下： H=2 -2;-2 4;f=-4;-12;A=1 1;-1 2;2 1;b=2;2;3;lb=zeros(2,1);x,fval,exitflag=quadprog(H,f,A,b,lb)三、非线性规划算法调用格式: x, fval, exitflag=fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlcon)说明：返回值 x为最优解向量，fval 是返回解所对应的目标函数值；exitflag 是描述搜索是否收敛。f 为目标函数，x0 为初
6、始点，A,b 为不等式约束的系数矩阵和右端列向量， 若没有不等式约束，则令 A= 、b= 。lb ,ub 为变量 x的下界和上界；nonlcon=fun，由 M文件 fun.m给定非线性不等式约束 c (x) 0 和非线性等式约束 g(x)=0。 例 5、求解非线性规划问题min 100(x2-x12 )2+(1-x1)2s.t x12;MATLAB 简介数学建模 x22首先建立 ff6.m文件：function f=ff6(x)f=100*(x(2)-x(1)2)2+(1-x(1)2;然后在命令窗口键入命令：x0=1.1,1.1; A=1 0;0 1;b=2;2;x,fval,exitfla
7、g=fmincon(ff6,x0,A,b)例 6、求解非线性规划问题min f=ex1(6x12+3x22+2x1x2+4x2+1)s.t x1x2-x1-x2+10-2x1x2-50首先建立目标函数文件 ff8.m文件：function f=ff8(x)f=exp(x(1)*(6*x(1)2+3*x(2)2+2*x(1)*x(2)+4*x(2)+1);再建立非线性的约束条件文件：ff8g.mfunction c,g=ff8g(x)c(1)=x(1)*x(2)-x(1)-x(2)+1;c(2)=-2*x(1)*x(2)-5;g=;然后在命令窗口键入以下命令：x0=1,1;nonlcon=ff8
8、gx,fval,exitflag =fmincon(ff8,x0, nonlcon)MATLAB 简介数学建模 四、整数线性规划算法说明：下面给出用分枝定界法求解整数线性规划的 M 函数文件 ILp.m，其中返回值 x为最优解向量，f 为最优值；x0 为初值点，可以用 代替；neqcstr 表示约束条件 Ax b中的前 neqcstr个是等式，neqcstr=0 时可以省略，此时也可以省略 x0；vlb ,vub 为变量 x的下界和上界；pre是精度。M 函数文件 ILp.m 如下：function x, f = ILp(c,A,b,vlb,vub,x0, neqcstr,pre)if nar
9、ginpre);mtemp=length(temp2);if isempty(temp2)x_f_b=xtemp; ftemp; vlb; vub;while jpre);if isempty(temp12),xall=xall,xtemp; fall=fall,ftemp;fvub=min(fvub,fall);elseif ftemppre);if isempty(tempr2),xall=xall,xtemp; fall=fall,ftemp;fvub=min(fvub,fall);elseif ftempmm,break,endtemp0=round (xint(:,j); temp1
10、=floor(xint(:,j);temp2=find(abs(xint(:,j)-temp0)pre);mtemp=length(temp2); end,else, x=xtemp; f=ftemp; end,if isempty(fall)fmin=min(fall); nmin=find(fall=fmin);x=xall(:,nmin); f=fmin; end,else, x=nan*ones(1,nvars); end例 7、求解整数线性规划问题max f=20x1+10x2s.t 5x1+4x2242x1+5x213xj0 ， i=1,2x1,x2 为整数先建立 M 函数文件 I
11、Lp.m，然后在 MATLAB 命令窗口键入：clear;c=-20,-10; %求 max 转换为求 mina=5,4;2,5; b=24;13;x,f=ILp(c,a,b,0;0,inf; inf, ,0,0.0001)f=-fMATLAB 简介数学建模 五、0-1 整数线性规划算法说明：下面的隐枚举法求解 01 线性规划的 M 函数文件 L01p_ie.m 中用到命令 B=de2bi(D)，其作用是将十进制数向量 D 转换为相应的二进制数按位构成的以 0，1 为元素的矩阵 B。M 函数文件 de2bi.m 如下：function b= de2bi(d,n,p)d=d(: ); len_d
12、=length(d);if min(d)0b1=b1 rem(tmp,2); tmp=floor(tmp/2);end;n=length(b1);end;if nargin0) b(i,j)=rem(tmp,p); tmp=floor(tmp/p);j=j+1;end;end;说明：下面给出用隐枚举法求解 01 线性规划的 M 函数文件 L01p_ie.m，其中返回值 x为最优解向量，f 为最优值；N 表示约束条件 Ax b中的前 N个是等式，N =0 时可以省略。M 函数文件 L01p_ie.m 如下：function x, f = L01p_ie(c,A,b,N)if nargin0, j
13、=1; end;endif j=m+1x=B; f=c*B;b(1)=min(b(1),f);endi=i+1;end例 8、求解 01整数线性规划问题max f=-3x1+2x2-5x3s.t x1+x2-x32x1+4x2+x34x1+x234x2+x36xj (j=1,2,3）为 0或 1先建立 M 函数文件 de2bi.m 和 L01p_ie.m，然后在 MATLAB 命令窗口键入：clear;c=3,-2,5; %求 max 转换为求 mina=1,2,-1;1,4,1;1,1,0;0,4,1; b=2;4;3;6;x, f = L01p_ie(c,a,b)f=-f以下整数线性规划算
14、法程序还未通过调试:M 函数文件 IntLP.m 如下：function x, y= IntLP(f,G,h,Geq,heq,lb,ub,x, id,options)global upper opt c x0 A b Aeq beq ID options;if nargin0.00005 %in order to avoid errorreturn;end;if max(abs(x.*ID-round(x.*ID)0.00005 %in order to avoid erroropt=x; upper=ftemp;return;elseopt=opt;x;return;end;end;noti
15、ntx=find(abs(x-round(x)=0.00005); %in order to avoid errorintx=fix(x); tempvlb=vlb; tempvub=vub;if vub(notintx(1,1),1)= intx(notintx(1,1),1)+1tempvlb(notintx(1,1),1)=intx(notintx(1,1),1)+1;ftemp=IntLP(tempvlb,vub);end;if vlb(notintx(1,1),1)= intx(notintx(1,1),1)tempvub(notintx(1,1),1)= intx(notintx(1,1),1);ftemp=IntLP(vlb, tempvub);end;clear;c=-20,-10; %求 max 转换为求 mina=5,4;2,5; b=24;13;x,y=IntLp(c,a,b, , ,0;0,inf; inf)f=-f
