多目标非线性规划程序Matlab完整版
- 格式:docx
- 大小:104.28 KB
- 文档页数:10
Matlab求解⾮线性规划,fmincon函数的⽤法总结Matlab求解⾮线性规划,fmincon函数的⽤法总结1.简介在matlab中,fmincon函数可以求解带约束的⾮线性多变量函数(Constrained nonlinear multivariable function)的最⼩值,即可以⽤来求解⾮线性规划问题matlab中,⾮线性规划模型的写法如下min\ f(x) \\ s.t. \begin{equation} \left\{ \begin{array}{**lr**} A \cdot x \leq b \\ Aeq\cdot x =beq\\ c(x)\leq0 \\ ceq(x)=0 \\ lb \leq x \leq ub\end{array} \right. \end{equation} \\ ~\\ f(x)是标量函数,x,b,beq是向量,A,Aeq是矩阵 \\ c(x)和ceq(x)是向量函数2.基本语法[x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)x的返回值是决策向量x的取值,fval的返回值是⽬标函数f(x)的取值fun是⽤M⽂件定义的函数f(x),代表了(⾮)线性⽬标函数x0是x的初始值A,b,Aeq,beq定义了线性约束 ,如果没有线性约束,则A=[],b=[],Aeq=[],beq=[]lb和ub是变量x的下界和上界,如果下界和上界没有约束,则lb=[],ub=[],也可以写成lb的各分量都为 -inf,ub的各分量都为infnonlcon是⽤M⽂件定义的⾮线性向量函数约束options定义了优化参数,不填写表⽰使⽤Matlab默认的参数设置3.实例⽰例,求下列⾮线性规划:min\ f(x)=x_1^2+x_2^2+x_3^2+8\\ s.t. \begin{equation} \left\{ \begin{array}{**lr**} x_1^2-x_2+x_3^2\geq0\\ x_1+x_2^2+x_3^2\leq20\\ -x_1-x_2^2+2=0\\ x_2+2x_3^2=3\\ x_1,x_2,x_3\geq0 \end{array} \right. \end{equation}(1)编写M函数fun1.m 定义⽬标函数:function f=fun1(x);f=x(1).^2+x(2).^2+x(3).^2+8;(2)编写M函数fun2.m定义⾮线性约束条件:function [g,h]=fun2(x);g=[-x(1).^2+x(2)-x(3).^2x(1)+x(2).^2+x(3).^3-20];h=[-x(1)-x(2).^2+2x(2)+2*x(3).^2-3];(3)编写主程序函数[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2')所得结果为:x_1=0.5522,x_2=1.2033,x_3=0.9478\\ 最⼩值y=10.651Processing math: 0%。
非线性整数规划的遗传算法Matlab程序(引自人工智能及数据挖掘论坛)这是一个具有200个01决策变量的多目标非线性整数规划,编写优化的目标函数如下,其中将多目标转化为单目标采用简单的加权处理。
function Fitness=FITNESS(x,FARM,e,q,w)%% 适应度函数% 输入参数列表% x 决策变量构成的4×50的0-1矩阵% FARM 细胞结构存储的当前种群,它包含了个体x% e 4×50的系数矩阵% q 4×50的系数矩阵% w 1×50的系数矩阵%%gamma=0.98;N=length(FARM);%种群规模F1=zeros(1,N);F2=zeros(1,N);for i=1:Nxx=FARM{i};ppp=(1-xx)+(1-q).*xx;F1(i)=sum(w.*prod(ppp));F2(i)=sum(sum(e.*xx));endppp=(1-x)+(1-q).*x;f1=sum(w.*prod(ppp));f2=sum(sum(e.*x));Fitness=gamma*sum(min([sign(f1-F1);zeros(1,N)]))+(1-gamma)*su m(min([sign(f2-F2);zeros(1,N)]));针对问题设计的遗传算法如下,其中对模型约束的处理是重点考虑的地方function [Xp,LC1,LC2,LC3,LC4]=MYGA(M,N,Pm)%% 求解01整数规划的遗传算法%% 输入参数列表% M 遗传进化迭代次数% N 种群规模% Pm 变异概率%% 输出参数列表% Xp 最优个体% LC1 子目标1的收敛曲线% LC2 子目标2的收敛曲线% LC3 平均适应度函数的收敛曲线% LC4 最优适应度函数的收敛曲线%% 参考调用格式[Xp,LC1,LC2,LC3,LC4]=MYGA(50,40,0.3) %% 第一步:载入数据和变量初始化load eqw;%载入三个系数矩阵e,q,w%输出变量初始化Xp=zeros(4,50);LC1=zeros(1,M);LC2=zeros(1,M);LC3=zeros(1,M);LC4=zeros(1,M);Best=inf;%% 第二步:随机产生初始种群farm=cell(1,N);%用于存储种群的细胞结构k=0;while k %以下是一个合法个体的产生过程x=zeros(4,50);%x每一列的1的个数随机决定for i=1:50R=rand;Col=zeros(4,1);if R<0.7RP=randperm(4);%1的位置也是随机的Col(RP(1))=1;elseif R>0.9RP=randperm(4);Col(RP(1:2))=1;elseRP=randperm(4);Col(RP(1:3))=1;endx(:,i)=Col;end%下面是检查行和是否满足约束的过程,对于不满足约束的予以抛弃Temp1=sum(x,2);Temp2=find(Temp1>20);if length(Temp2)==0k=k+1;farm{k}=x;endend%% 以下是进化迭代过程counter=0;%设置迭代计数器while counter%% 第三步:交叉%交叉采用双亲双子单点交叉newfarm=cell(1,2*N);%用于存储子代的细胞结构Ser=randperm(N);%两两随机配对的配对表A=farm{Ser(1)};%取出父代AB=farm{Ser(2)};%取出父代BP0=unidrnd(49);%随机选择交叉点a=[A(:,10),B(:,(P0+1):end)];%产生子代ab=[B(:,10),A(:,(P0+1):end)];%产生子代bnewfarm{2*N-1}=a;%加入子代种群newfarm{2*N}=b;%以下循环是重复上述过程for i=1N-1)A=farm{Ser(i)};B=farm{Ser(i+1)};P0=unidrnd(49);a=[A(:,10),B(:,(P0+1):end)];b=[B(:,10),A(:,(P0+1):end)];newfarm{2*i-1}=a;newfarm{2*i}=b;endFARM=[farm,newfarm];%新旧种群合并%% 第四步:选择复制FLAG=ones(1,3*N);%标志向量,对是否满足约束进行标记%以下过程是检测新个体是否满足约束for i=13*N)x=FARM{i};sum1=sum(x,1);sum2=sum(x,2);flag1=find(sum1==0);flag2=find(sum1==4);flag3=find(sum2>20);if length(flag1)+length(flag2)+length(flag3)>0FLAG(i)=0;%如果不满足约束,用0加以标记endendNN=length(find(FLAG)==1);%满足约束的个体数目,它一定大于等于N NEWFARM=cell(1,NN);%以下过程是剔除不满主约束的个体kk=0;for i=13*N)。
多元⾮线性规划Matlab,⾮线性规划MATLAB代码下⾯是三个⾮线性规划领域的算法。
课堂上给予了详细的讲解,在实践环节让学⽣编程实现,从⽽可以实验复杂⼀些的例⼦,加深对算法的理解。
下⾯共有四个程序grad,simplelinesearch,bfgs和phr,全部使⽤MATLAB语⾔编写。
这些代码远未完善,可修改余地很⼤,仅供教学之⽤。
function gradf=grad(hfun,x)%GRAD 数值法求函数在给定点处的导数值(⼀元函数)或梯度(多元函数)% gradf = grad(hfun,x0) hfun是函数句柄或内联函数,x0是⼀定点或⼀批点(按列);返回% 值gradf是函数在该点处的导数或梯度。
% 要求函数能对成批的点求函数值。
⽐如:feval(hfun,X)返回⼀个与X同列的⾏向量,对应于以% X每⼀列作为函数⾃变量⽽求得的函数值。
%% Reference: 《最优化计算原理与算法程序设计》, 粟塔⼭等编著, 国防科技⼤学出版社% (湖南), 2001.%% $Author: WBC $ $Date: 2003/10/25 $n = length(x);h=1e-3; % 数值法求梯度的步长w1=zeros(n,1);%h/2w2=w1; %-h/2w3=w1; %hw4=w1; %-hfor i=1:nx(i)=x(i)+h/2;w1(i)=hfun(x);x(i)=x(i)-h;w2(i)=hfun(x);x(i)=x(i)-h/2;w4(i)=hfun(x);x(i)=x(i)+2*h;w3(i)=hfun(x);x(i)=x(i)-h;endgradf=(8*(w1-w2)-(w3-w4))/(6*h);function [fv,x,lambda,exitflag]=simplelinesearch(hf,rho,l,u,lambda0,fv0,x0,g,d)%LINESEARCH 简单线搜索%输⼊参数:% hf -- 函数句柄,⽬标函数% rho -- 实标量,简单线搜索的参数,介于0到0.5之间的数% l -- 实标量,简单线搜索的参数,介于0到1之间的数% u -- 实标量,简单线搜索的参数,介于0到1之间的数,满⾜u>l% alpha0 -- 实标量,步长的初始值% fv0 -- 实标量,⼀维搜索的初始⽬标函数值,即 hf(x0+alpha_0*d)% x0 -- 实列向量,当前点% g -- 实列向量,函数hf在当前点处的梯度% d -- 实列向量,函数hf在当前点处的搜索⽅向%输出参数:% fv -- 实标量,⼀维搜索完成后的⽬标函数值,即 hf(x0+alpha*d)% x -- 实列向量,下⼀个点% lambda -- 实标量,可接受步长% exitflag -- 整型标量,等于0表⽰线搜索成功,等于-1表⽰线搜索失败(内部迭代次数⼤于iterMax) %参考⽂献:倪勤,最优化⽅法与程序设计,科学出版社% Date: 2009/12/20lambda = lambda0;x = x0;i = 0;imax = 30; % 最⼤迭代次数,⽤户可以修改% 主循环gd = dot(g, d);while i <= imaxfv = hf(x + lambda*d);i = i + 1;if fv < fv0 + lambda * rho * gd;x = x + lambda*d;exitflag = 0;return;endlambda_bar = -gd*lambda^2*0.5/(fv-fv0-lambda*gd);lambda = min(lambda_bar, u*lambda);endexitflag=0;if i >= imax && fv >= fv0fv = fv0;x = x0;lambda = 0;exitflag = -1;endfunction [fv, x, exitflag] = bfgs(hf, x0, epsi)%BFGS ⽆约束问题的BFGS算法%输⼊参数:% hf -- 函数句柄,⽬标函数% x0 -- 实列向量,初始点% epsi -- 实标量,终⽌误差%输出参数:% fv -- 实标量,⼀维搜索完成后的⽬标函数值,即 hf(x0+alpha*d)% x -- 实列向量,下⼀个点% exitflag -- 整型标量,等于0表⽰成功,等于-1表⽰失败(迭代次数⼤于iter_max) %参考⽂献:倪勤,最优化⽅法与程序设计,科学出版社% Date: 2009/12/20%%初始化%k = 0;rho = 0.01;l = 0.15;u = 0.85;x =x0;fv = hf(x);n = length(x);H = eye(n);iter_max = 100;%检查终⽌条件%g = grad(hf, x);while norm(g) > epsi && k<= iter_maxd = -H*g;lambda0 = 1.0;%%做线搜索,如果成功,则返回更新当前点%[fv, x, lambda, exitflag] = simplelinesearch(hf, rho, l, u, lambda0, fv, x, g, d); if exitflag == -1 %重开始H = eye(n);else%%更新H%g_old = g;g = grad(hf, x);p = lambda*d;q = g- g_old;Hq = H*q;pq = p'*q;qHq = q'*Hq;v = sqrt(qHq) * (p/pq-Hq/qHq);H = H + p*p'/pq - Hq*Hq'/qHq + v*v';endk = k+1;endexitflag = 0;if k > iter_maxexitflag = -1;endfunction [fv, x, exitflag] = phr(hf, cf, x0)%输⼊参数:% hf -- 函数句柄,⽬标函数% cf -- 函数句柄,约束条件,包含等式约束和不等式约束% x0 -- 实列向量,初始点%输出参数:% fv -- 实标量,⼀维搜索完成后的⽬标函数值,即 hf(x0+alpha*d)% x -- 实列向量,下⼀个点% exitflag -- 整型标量,等于0表⽰成功,等于-1表⽰失败(迭代次数⼤于iter_max) %参考⽂献:倪勤,最优化⽅法与程序设计,科学出版社% Date: 2009/12/20%%初始化%epsi = 1.0e-4;k = 0;sigma = 0.8;c = 1.5;theta = 0.8;x = x0;[ce, ci] = cf(x);l = length(ce);li = length(ci);lambda = ones(l+li, 1) * 0.1;iter_max = 100;phi = 0;if lphi = phi + ce'*ce;endif liphi = phi + sum(min(ci,lambda(l+1:end)/sigma).^2);endwhile phi > epsi && k <= iter_max%%hmf = @(x) mfun(x, hf, cf,lambda, sigma);[fv, x] = bfgs(hmf, x, epsi);[ce, ci] = cf(x);phi_old = phi;phi = 0;if lphi = phi + ce'*ce;endif liphi = phi + sum(min(ci,lambda(l+1:end)/sigma).^2); endif phi > epsi%%更新罚因⼦%if k >= 2 && phi/phi_old > thetasigma = c * sigma;end%%更新乘⼦%if llambda(1:l) = lambda(1:l) - sigma*ce;endif lilambda(l+1:end) = max(0, lambda(l+1:end) - sigma*ci); endendk = k+1;endexitflag = 0;fv = hf(x);exitflag = -1;end%%乘⼦罚函数%function fv = mfun(x, hf, cf, lambda, sigma)[ce, ci] = cf(x);l = length(ce);li = length(ci);fv = 0;fv = fv + hf(x);if lfv = fv - lambda(1:l)'*ce + 0.5*sigma*ce'*ce;endif lifv = fv + 0.5/sigma*sum(max(0,lambda(l+1:end) - sigma*ci).^2 - lambda(l+1:end).^2); end这⾥是演⽰代码:%% bfgs演⽰%教材P328.1-3hf1 = @(x) 100 * (x(2) - x(1).^2).^2 + (1 - x(1)).^2; %banana函数hf2 = @(x) (6 + x(1) + x(2)).^2 + (2 - 3*x(1) - 3*x(2) - x(1)*x(2)).^2;hf3 = @(x) x(1).^2 - 2*x(1)*x(2) + 4*x(2).^2 + x(1) - 3*x(2);[fv,x,exitflag]=bfgs(hf1,[0;0],0.001);fvxexitflag[fv,x,exitflag]=bfgs(hf2,[4;6],0.001);fvxexitflag[fv,x,exitflag]=bfgs(hf3,[1;1],0.001);fvexitflag%% phr算法%教材P414.5hf1 = @(x) x(1).^2 + x(2).^2;cf1 = @(x) deal([], x(1) - 1);hf2 = @(x) x(1) + (x(2) + 1).^2/3;cf2 = @(x) deal([],[x(1); x(2) - 1]);%教材P392.2hf3 = @(x) x(1).^2 + x(1).*x(2) + 2*x(2).^2 - 6*x(1) - 2*x(2) - 12*x(3);cf3 = @(x) deal(x(1)+x(2)+x(3)-2,...[x(1) - 2*x(2) + 3; x(1); x(2); x(3)]);%其它例⼦hf4 = @(x) 6*x(2)*x(5) + 7*x(1)*x(3) + 3*x(2)^2;cf4 = @(x) deal([3*x(2)^2*x(5) + 3*x(1)^2*x(3) - 20.875;x(1) - 0.3*x(2)],...[-x(1) + 0.2*x(2)*x(5) + 71-0.9*x(3) + x(4)^2 + 67x(3)x(5) - 1-x(3) + 20x(4) - 0.1*x(5)-x(4) + 0.5*x(5)x(3) - 0.9*x(5)]);hf5 = @(x) exp(x(1)) * (4*x(1)^2 + 2*x(2)^2 + 4*x(1)*x(2) + 2*x(2) + 1);cf5 = @(x) deal([], [x(1) + x(2) - x(1)*x(2) - 1.5; x(1)*x(2) + 10]);[fv, x, exitflag] = phr(hf1, cf1, [3;2]); fv x exitflag [fv, x, exitflag] = phr(hf2, cf2, [3;2]); fv x exitflag [fv, x, exitflag] =phr(hf3, cf3, [1;1;0]); fv x exitflag [fv, x, exitflag] = phr(hf4, cf4, [1; 4; 5; 2; 5]); fv x exitflag [fv, x, exitflag] = phr(hf5, cf5, [-1; 1]); fv x exitflag。
非线性整数规划的遗传算法Matlab程序(附图)通常,非线性整数规划是一个具有指数复杂度的NP问题,如果约束较为复杂,Matlab 优化工具箱和一些优化软件比如lingo等,常常无法应用,即使能应用也不能给出一个较为令人满意的解。
这时就需要针对问题设计专门的优化算法。
下面举一个遗传算法应用于非线性整数规划的编程实例,供大家参考!模型的形式和适应度函数定义如下:这是一个具有200个01决策变量的多目标非线性整数规划,编写优化的目标函数如下,其中将多目标转化为单目标采用简单的加权处理。
function Fitness=FITNESS(x,FARM,e,q,w)%% 适应度函数% 输入参数列表% x 决策变量构成的4×50的0-1矩阵% FARM 细胞结构存储的当前种群,它包含了个体x% e 4×50的系数矩阵% q 4×50的系数矩阵% w 1×50的系数矩阵%%gamma=0.98;N=length(FARM);%种群规模F1=zeros(1,N);F2=zeros(1,N);for i=1:Nxx=FARM{i};ppp=(1-xx)+(1-q).*xx;F1(i)=sum(w.*prod(ppp));F2(i)=sum(sum(e.*xx));endppp=(1-x)+(1-q).*x;f1=sum(w.*prod(ppp));f2=sum(sum(e.*x));Fitness=gamma*sum(min([sign(f1-F1);zeros(1,N)]))+(1-gamma)*sum(mi n([sign(f2-F2);zeros(1,N)]));针对问题设计的遗传算法如下,其中对模型约束的处理是重点考虑的地方function [Xp,LC1,LC2,LC3,LC4]=MYGA(M,N,Pm)%% 求解01整数规划的遗传算法%% 输入参数列表% M 遗传进化迭代次数% N 种群规模% Pm 变异概率%% 输出参数列表% Xp 最优个体% LC1 子目标1的收敛曲线% LC2 子目标2的收敛曲线% LC3 平均适应度函数的收敛曲线% LC4 最优适应度函数的收敛曲线%% 参考调用格式[Xp,LC1,LC2,LC3,LC4]=MYGA(50,40,0.3)%% 第一步:载入数据和变量初始化load eqw;%载入三个系数矩阵e,q,w%输出变量初始化Xp=zeros(4,50);LC1=zeros(1,M);LC2=zeros(1,M);LC3=zeros(1,M);LC4=zeros(1,M);Best=inf;%% 第二步:随机产生初始种群farm=cell(1,N);%用于存储种群的细胞结构k=0;while k %以下是一个合法个体的产生过程x=zeros(4,50);%x每一列的1的个数随机决定for i=1:50R=rand;Col=zeros(4,1);if R<0.7RP=randperm(4);%1的位置也是随机的Col(RP(1))=1;elseif R>0.9RP=randperm(4);Col(RP(1:2))=1;elseRP=randperm(4);Col(RP(1:3))=1;endx(:,i)=Col;end%下面是检查行和是否满足约束的过程,对于不满足约束的予以抛弃 Temp1=sum(x,2);Temp2=find(Temp1>20);if length(Temp2)==0k=k+1;farm{k}=x;endend%% 以下是进化迭代过程counter=0;%设置迭代计数器while counter%% 第三步:交叉%交叉采用双亲双子单点交叉newfarm=cell(1,2*N);%用于存储子代的细胞结构Ser=randperm(N);%两两随机配对的配对表A=farm{Ser(1)};%取出父代AB=farm{Ser(2)};%取出父代BP0=unidrnd(49);%随机选择交叉点a=[A(:,1:P0),B(:,(P0+1):end)];%产生子代ab=[B(:,1:P0),A(:,(P0+1):end)];%产生子代bnewfarm{2*N-1}=a;%加入子代种群newfarm{2*N}=b;%以下循环是重复上述过程for i=1:(N-1)A=farm{Ser(i)};B=farm{Ser(i+1)};P0=unidrnd(49);a=[A(:,1:P0),B(:,(P0+1):end)];b=[B(:,1:P0),A(:,(P0+1):end)];newfarm{2*i-1}=a;newfarm{2*i}=b;endFARM=[farm,newfarm];%新旧种群合并%% 第四步:选择复制FLAG=ones(1,3*N);%标志向量,对是否满足约束进行标记%以下过程是检测新个体是否满足约束for i=1:(3*N)x=FARM{i};sum1=sum(x,1);sum2=sum(x,2);flag1=find(sum1==0);flag2=find(sum1==4);flag3=find(sum2>20);if length(flag1)+length(flag2)+length(flag3)>0FLAG(i)=0;%如果不满足约束,用0加以标记endendNN=length(find(FLAG)==1);%满足约束的个体数目,它一定大于等于N NEWFARM=cell(1,NN);%以下过程是剔除不满主约束的个体kk=0;for i=1:(3*N)if FLAG(i)==1kk=kk+1;NEWFARM{kk}=FARM{i};endend%以下过程是计算并存储当前种群每个个体的适应值SYZ=zeros(1,NN);syz=zeros(1,N);for i=1:NNx=NEWFARM{i};SYZ(i)=FITNESS2(x,NEWFARM,e,q,w);%调用适应值子函数endk=0;%下面是选择复制,选择较优的N个个体复制到下一代while k minSYZ=min(SYZ);posSYZ=find(SYZ==minSYZ);POS=posSYZ(1);k=k+1;farm{k}=NEWFARM{POS};syz(k)=SYZ(POS);SYZ(POS)=inf;end%记录和更新,更新最优个体,记录收敛曲线的数据minsyz=min(syz);meansyz=mean(syz);pos=find(syz==minsyz);LC3(counter+1)=meansyz;if minsyz Best=minsyz;Xp=farm{pos(1)};endLC4(counter+1)=Best;ppp=(1-Xp)+(1-q).*Xp;LC1(counter+1)=sum(w.*prod(ppp));LC2(counter+1)=sum(sum(e.*Xp));%% 第五步:变异for i=1:Nif Pm>rand%是否变异由变异概率Pm控制AA=farm{i};%取出一个个体POS=unidrnd(50);%随机选择变异位R=rand;Col=zeros(4,1);if R<0.7RP=randperm(4);Col(RP(1))=1;elseif R>0.9RP=randperm(4);Col(RP(1:2))=1;elseRP=randperm(4);Col(RP(1:3))=1;end%下面是判断变异产生的新个体是否满足约束,如果不满足,此次变异无效 AA(:,POS)=Col;Temp1=sum(AA,2);Temp2=find(Temp1>20);if length(Temp2)==0farm{i}=AA;endendendcounter=counter+1end%第七步:绘收敛曲线图figure(1);plot(LC1);xlabel('迭代次数');ylabel('子目标1的值');title('子目标1的收敛曲线'); figure(2);plot(LC2);xlabel('迭代次数');ylabel('子目标2的值');title('子目标2的收敛曲线'); figure(3);plot(LC3);xlabel('迭代次数');ylabel('适应度函数的平均值');title('平均适应度函数的收敛曲线'); figure(4);plot(LC4);xlabel('迭代次数');ylabel('适应度函数的最优值');title('最优适应度函数的收敛曲线');贴出一幅运行得到的收敛曲线。
在matlab 中非线性规划的数学模型可写成一下形式:minf(X)s.t. Ax ≪B Aeq .x =Beq C (x )≪0Ceq x =0其中,f(x)是标量函数;A,B,Aeq,Beq 是相应维数的矩阵和向量;C(x),Ceq(x)是非线性向量函数。
Matlab 中的命令是X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)它的返回值是向量x 。
其中,FUN 是用M 文件定义的函数f(x)。
X0是X 的初始值。
A ,B ,Aeq ,Beq 定义了线性约束AX ≪B ,Aeq*X=Beq ,如果没有线性约束,则A=[],B=[],Aeq=[],Beq=[]。
LB 和UB 是变量x 的下界和上界,如果上界和下界没有约束,则LB=[],UB=[];如果X 无下界,则LB=-inf;如果X 无上界,则UB=inf 。
NONLCON 是用M 文件定义的非线性向量函数C(x),Ceq(x)。
OPTIONS 定义了优化函数,可以使用MATLAB 默认的参数设置。
例求解下列非线性规划问题:max z= X 1+ X 2+ X 3+ X 4 s.t.x 1≪4001.1x 1+x 2≪4401.21x 1+1.1x 2+x 3≪4841.331x 1+1.21x 2+1.1x 3+x 4≪532.4X i≫0,i =1,2,3,4(1)编写M 文件,定义目标函数:function f=fun44(x)f=-(sqrt(x(1))+sqrt(x(2))+sqrt(x(3))+sqrt(x(4)) );(2)编写M 文件,定义约束条件function[g,ceq]=mycon1(x)g(1)=x(1)-400;g(2)=1.1*x(1)+x(2)-440;g(3)=1.21*x(1)+1.1*x(2)+x(3)-484;g(4)=1.331*x(1)+1.21*x(2)+1.1*x(3)+x(4)-532.4;ceq=0(3)编写主程序x0=[1;1;1;1];lb=[0;0;0;0];ub=[];A=[];b=[];Aeq=[];beq=[];[x,fval] = fmincon('fun44',x0,A,b,Aeq,beq,lb,ub,'mycon1')输出结果x =86.1883104.2879 126.1883 152.6879fval =-43.0860。
Matlab 数学建模学习笔记——⾮线性规划⽬录⾮线性规划⾮线性规划的matlab 解法fmincon 函数x = fmincon(fun,x0,A,b)x = fmincon(fun,x0,A,b,Aeq,beq)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)% x0是x 的初始值,fun 是⽤M ⽂件定义的函数f(x),nonlcon 是M ⽂件定义的⾮线性向量函数c(x),ceq(x);options 定义了优化参数,可以⽤Matlab 默认的参数设置x = fmincon(problem)[x,fval] = fmincon(___)[x,fval,exitflag,output] = fmincon(___)[x,fval,exitflag,output,lambda,grad,hessian] = fmincon(___)e.gmin f (x )=x 21+x 22+x 23+8,s .t .=x 21−x 2+x 23≥0,x 1+x 22+x 33≤20,−x 1−x 22+2=0,x 2+2x 23=3,x 1,x 2,x 3≥0。
对于这道题,我们需要编写三个m ⽂件% filename :'fun1.m'function f = fun1(x)f = sum(x.^2)+8;% filename :'fun2.m'function [g,h] = fun2(x)g = [-x(1)^2+x(2)-x(3)^2x(1)+x(2)^2+x(3)^3-20]; % ⾮线性不等式约束h = [-x(1)-x(2)^2+2x(2)+2*x(3)^2-3]; % ⾮线性等式约束% filename :'main.m'clear;clc;[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[],'fun2');求得当x1=0.5522, x2=1.2033, x3=0.9478时,最⼩值y=10.6511fminsearch 函数(求极⼩值)f (x )=sin (x )+3例如我们要求f (x )取得极⼩值时x 的值,代码如下% filename :'fun.m'function f = fun(x);f = sin(x) +3;% filename :'main.m'x0 = 2;[x,y] = fminsearch(@ fun,x0)计算函数的零点和⽅程组的解求多项式f (x )=x 3−x 2+2x −3的零点法⼀clear;clc;xishu = [1,-1,2,-3];x = roots(xishu);x =-0.1378 + 1.5273i -0.1378 - 1.5273i 1.2757 + 0.0000i 法⼆clear;clc;syms xx0 = solve(x^3-x^2+2*x-3);%求函数零点的符号解x0 = vpa(x0,5); % 化成⼩数格式的数据,5为有效数字法三clear;clc;y = @(x)x^3-x^2+2*x-3;x = fsolve(y,rand);% 只能求给定初始值附近的⼀个零点约束极值问题{Processing math: 100%⼆次规划H为实对称矩阵求解命令x = quadprog(H,f)x = quadprog(H,f,A,b)x = quadprog(H,f,A,b,Aeq,beq)x = quadprog(H,f,A,b,Aeq,beq,lb,ub)x = quadprog(H,f,A,b,Aeq,beq,lb,ub,x0)x = quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,options)x = quadprog(problem)[x,fval] = quadprog(___)[x,fval,exitflag,output] = quadprog(___)[x,fval,exitflag,output,lambda] = quadprog(___)罚函数法利⽤罚函数法,可将⾮线性规划规划问题的求解,转化为求解⼀系列⽆约束极值问题,也可称这种⽅法为序列⽆约束最⼩化技术思想是利⽤题⽬中的约束函数做出适当的罚函数,由此构造出带参数的增⼴⽬标函数,把问题转化成⽆约束⾮线性规划问题。
MATLAB⾮线性规划问题⼀.⾮线性规划课题实例1 表⾯积为36平⽅⽶的最⼤长⽅体体积。
建⽴数学模型:设x、y、z分别为长⽅体的三个棱长,f为长⽅体体积。
max f = x y (36-2 x y)/2 (x+y)实例2 投资决策问题某公司准备⽤5000万元⽤于A、B两个项⽬的投资,设x1、x2分别表⽰配给项⽬A、B的投资。
预计项⽬A、B的年收益分别为20%和16%。
同时,投资后总的风险损失将随着总投资和单位投资的增加⽽增加,已知总的风险损失为2x12+x22+(x1+x2)2.问应如何分配资⾦,才能使期望的收益最⼤,同时使风险损失为最⼩。
建⽴数学模型:max f=20x1+16x2-λ[2x12+x22+(x1+x2)2]s.t x1+x2≤5000x 1≥0,x2≥0⽬标函数中的λ≥0是权重系数。
由以上实例去掉实际背景,其⽬标函数与约束条件⾄少有⼀处是⾮线性的,称其为⾮线性问题。
⾮线性规划问题可分为⽆约束问题和有约束问题。
实例1为⽆约束问题,实例2为有约束问题。
⼆.⽆约束⾮线性规划问题:求解⽆约束最优化问题的⽅法主要有两类:直接搜索法(Search method)和梯度法(Gradient method),单变量⽤fminbnd,fminsearch,fminunc;多变量⽤fminsearch,fminnuc 1.fminunc函数调⽤格式:x=fminunc(fun,x0)x=fminunc(fun,x0,options)x=fminunc(fun,x0,options,P1,P2)[x,fval]=fminunc(…)[x,fval, exitflag]=fminunc(…)[x,fval, exitflag,output]=fminunc(…)[x,fval, exitflag,output,grad]=fminunc(…)[x,fval, exitflag,output,grad,hessian]=fminunc(…)说明:fun为需最⼩化的⽬标函数,x0为给定的搜索的初始点。
多目标非线性规划程序M a t l a bDocument serial number【NL89WT-NY98YT-NC8CB-NNUUT-NUT108】f u n c t i o n[e r r m s g,Z,X,t,c,f a i l]=BNB18(fun,x0,xstat,xl,xu,A,B,Aeq,Beq,nonlcon,setts,options1,options2,maxSQPit,varargin );%·Dêy1£Díóa·§¨μü′ú·¨£úDê1ó£DèOptimization toolbox §3% Minimize F(x)%subject to: xlb <= x <=xub% A*x <= B% Aeq*x=Beq% C(x)<=0% Ceq(x)=0%% x(i)éaáD±á£êy£ò1ì¨μ% ê1óê%[errmsg,Z,X]=BNB18('fun',x0,xstat,xl,xu,A,B,Aeq,Beq,'nonlcon',setts)%fun£o Mt£±íê×Dˉ±êoˉêyf=fun(x)%x0: áDòᣱíê±á3μ%xstat£o áDòá£xstat(i)=0±íêx(i)aáD±á£1±íêêy£2±íê1ì¨μ%xl£o áDòᣱíê±á%xu: áDòᣱíê±áé%A: ó, ±íêD2μèêêμêy%B: áDòá, ±íêD2μèêêé%Aeq: ó, ±íêDμèêêμêy%Beg: áDòá, ±íêD2μèêêóòμ%nonlcon: Mt£±íê·Dêoˉêy[C,Ceq]=nonlin(x),DC(x)a2μèêê,% Ceq(x)aμèêê%setts: ·¨éè%errmsq: ·μ′íóìáê%Z: ·μ±êoˉêy×Dμ%X: ·μ×óa%%àyìa% max x1*x2*x3% -x1+2*x2+2*x3>=0% x1+2*x2+2*x3<=72% 10<=x2<=20% x1-x2=10% èD′ Moˉêy% function f=discfun(x)% f=-x(1)*x(2)*x(3);%óa% clear;x0=[25,15,10]';xstat=[1 1 1]';% xl=[20 10 -10]';xu=[30 20 20]';% A=[1 -2 -2;1 2 2];B=[0 72]';Aeq=[1 -1 0];Beq=10;% [err,Z,X]=BNB18('discfun',x0,xstat,xl,xu,A,B,Aeq,Beq);% XMAX=X',ZMAX=-Z%% BNB18 Finds the constrained minimum of a function of several possibly integer variables.% Usage: [errmsg,Z,X,t,c,fail] =%BNB18(fun,x0,xstatus,xlb,xub,A,B,Aeq,Beq,nonlcon,settings,options1,options2,maxSQPiter ,P1,P2,...)%% BNB solves problems of the form:% Minimize F(x) subject to: xlb <= x0 <=xub% A*x <= B Aeq*x=Beq% C(x)<=0 Ceq(x)=0% x(i) is continuous for xstatus(i)=0% x(i) integer for xstatus(i)= 1% x(i) fixed for xstatus(i)=2%% BNB uses:% Optimization Toolbox Version (R11) 09-Oct-1998% From this toolbox is called. For more info type help fmincon.%% fun is the function to be minimized and should return a scalar. F(x)=feval(fun,x).% x0 is the starting point for x. x0 should be a column vector.% xstatus is a column vector describing the status of every variable x(i).% xlb and xub are column vectors with lower and upper bounds for x.% A and Aeq are matrices for the linear constrains.% B and Beq are column vectors for the linear constrains.% nonlcon is the function for the nonlinear constrains.% [C(x);Ceq(x)]=feval(nonlcon,x). Both C(x) and Ceq(x) should be column vectors.%% errmsg is a string containing an error message if BNB found an error in the input.% Z is the scalar result of the minimization, X the values of the accompanying variables.% t is the time elapsed while the algorithm BNB has run, c is the number of BNB cycles and% fail is the number of unsolved leaf sub-problems.%% settings is a row vector with settings for BNB:% settings(1) (standard 0) if 1: use phase 1 by relaxation. This sometimes makes the algorithm% faster, because phase 1 means the algorithm first checks if there is a feasible solution% for a sub-problem before trying to find a best solution. If there is no feasible solution BNB% will not try to find a best solution.% settings(2) (standard 0) if 1: if the sub-problem did not converge do not branch. If a sub-% problem did not converge this means BNB did not find a solution for it. Normally BNB will% branch the problem so it can try again to find a solution.% A sub-problem that is a leaf of the branch-and-bound-three can not be branched. If such% a problem does not converge it will be considered unfeasible and the parameter fail will be% raised by one.% settings(3) (standard 0) if 1: if 1 a sub-problem that did not converge but did return a feasible% point will be considered convergent. This might be useful if fmincon is having a hard time with% a certain problem but you do want some results.% options1 and options2 are options structures for phase 1 and phase 2.% For details about the options structure type help optimset.% maxSQPiter is a global variable used by fmincon (if modified as described in .% maxSQPiter is 1000 by default.% P1,P2,... are parameters to be passed to fun and nonlcon.% F(x)=feval(fun,x,P1,P2,...). [C(x);Ceq(x)]=feval(nonlcon,x,P1,P2,...).% Type edit BNB18 for more info.% . Kuipers% e-mail% FI-Lab% Applied Physics% Rijksuniversiteit Groningen% To get rid of bugs and to stop fmincon from hanging make the following chances:%% In optim/private/ ($Revision: $ $Date: 1998/08/24 13:46:15 $):% Get EXITFLAG independent of verbosity.% After the lines: disp(' less than 2* but constraints are not satisfied.')% end% EXITFLAG = -1;% end% end% status=1;% add the line: if (strncmp(howqp, 'i',1) & mg > 0), EXITFLAG = -1; end;%% In optim/private/ ($Revision: $ $Date: 1998/09/01 21:37:56 $):% Stop qpsub from hanging.% After the line: % Andy Grace 7-9-90. Mary Ann Branch 9-30-96.% add the line: global maxSQPiter;% and changed the line: maxSQPiters = Inf;% to the line: if exist('maxSQPiter','var'), maxSQPiters = maxSQPiter; else maxSQPiters=inf; end;% I guess there was a reason to put maxSQPiters at infinity, but this works fine for me.global maxSQPiter;% STEP 0 CHECKING INPUTZ=[]; X=[]; t=0; c=0; fail=0;if nargin<2, errmsg='BNB needs at least 2 input arguments.'; return; end;if isempty(fun), errmsg='No fun found.'; return; end;if isempty(x0), errmsg='No x0 found.'; return;elseif size(x0,2)>1, errmsg='x0 must be a column vector.'; return; end;xstatus=zeros(size(x0));if nargin>2 & ~isempty(xstat)if all(size(xstat)<=size(x0))xstatus(1:size(xstat))=xstat;else errmsg='xstatus must be a column vector the same size as x0.'; return;end;if any(xstatus~=round(xstatus) | xstatus<0 | 2<xstatus)errmsg='xstatus must consist of the integers 0,1 en 2.'; return;end;end;xlb=zeros(size(x0));xlb(find(xstatus==0))=-inf;if nargin>3 & ~isempty(xl)if all(size(xl)<=size(x0))xlb(1:size(xl,1))=xl;else errmsg='xlb must be a column vector the same size as x0.'; return;end;end;if any(x0<xlb)errmsg='x0 must be in the range xlb <= x0.'; return;elseif any(xstatus==1 & (~isfinite(xlb) | xlb~=round(xlb)))errmsg='xlb(i) must be an integer if x(i) is an integer variabele.'; return;end;xlb(find(xstatus==2))=x0(find(xstatus==2));xub=ones(size(x0));xub(find(xstatus==0))=inf;if nargin>4 & ~isempty(xu)if all(size(xu)<=size(x0))xub(1:size(xu,1))=xu;else errmsg='xub must be a column vector the same size as x0.'; return;end;end;if any(x0>xub)errmsg='x0 must be in the range x0 <=xub.'; return;elseif any(xstatus==1 & (~isfinite(xub) | xub~=round(xub)))errmsg='xub(i) must be an integer if x(i) is an integer variabale.'; return;end;xub(find(xstatus==2))=x0(find(xstatus==2));if nargin>5if ~isempty(A) & size(A,2)~=size(x0,1), errmsg='Matrix A not correct.'; return; end; else A=[]; end;if nargin>6if ~isempty(B) & any(size(B)~=[size(A,1) 1]), errmsg='Column vector B not correct.'; return; end;else B=[]; end;if isempty(A) & ~isempty(B), errmsg='A and B should only be nonempty together.'; return; end;if isempty(B) & ~isempty(A), B=zeros(size(A,1),1); end;if nargin>7 & ~isempty(Aeq)if size(Aeq,2)~=size(x0,1), errmsg='Matrix Aeq not correct.'; return; end;else Aeq=[]; end;if nargin>8if ~isempty(Beq) & any(size(Beq)~=[size(Aeq,1) 1]), errmsg='Column vector Beq not correct.'; return; end;else Beq=[]; end;if isempty(Aeq) & ~isempty(Beq), errmsg='Aeq and Beq should only be nonempty together'; return; end;if isempty(Beq) & ~isempty(Aeq), Beq=zeros(size(Aeq,1),1); end;if nargin<10, nonlcon=''; end;settings = [0 0 0];if nargin>10 & ~isempty(setts)if all(size(setts)<=size(settings))settings(setts~=0)=setts(setts~=0);else errmsg='settings should be a row vector of length 3.'; return; end;end;if nargin<12, options1=[]; end;options1=optimset(optimset('fmincon'),options1);if nargin<13, options2=[]; end;options2=optimset(optimset('fmincon'),options2);if nargin<14, maxSQPiter=1000;elseif isnumeric(maxSQPit) & all(size(maxSQPit))==1 & maxSQPit>0 &round(maxSQPit)==maxSQPitmaxSQPiter=maxSQPit;else errmsg='maxSQPiter must be an integer >0'; return; end;eval(['z=',fun,'(x0,varargin{:});'],'errmsg=''fun caused error.''; return;');if ~isempty(nonlcon)eval(['[C, Ceq]=',nonlcon,'(x0,varargin{:});'],'errmsg=''nonlcon caused error.''; return;');if size(C,2)>1 | size(Ceq,2)>1, errmsg='C en Ceq must be column vectors.'; return; end;end;% STEP 1 INITIALISATIONcurrentwarningstate=warning;warning off;tic;lx = size(x0,1);z_incumbent=inf;x_incumbent=inf*ones(size(x0));I = ceil(sum(log2(xub(find(xstatus==1))-xlb(find(xstatus==1))+1))+size(find(xstatus==1),1)+1);stackx0=zeros(lx,I);stackx0(:,1)=x0;stackxlb=zeros(lx,I);stackxlb(:,1)=xlb;stackxub=zeros(lx,I);stackxub(:,1)=xub;stacksize=1;xchoice=zeros(size(x0));if ~isempty(Aeq)j=0;for i=1:size(Aeq,1)if Beq(i)==1 & all(Aeq(i,:)==0 | Aeq(i,:)==1)J=find(Aeq(i,:)==1);if all(xstatus(J)~=0 & xchoice(J)==0 & xlb(J)==0 & xub(J)==1) if all(xstatus(J)~=2) | all(x0(J(find(xstatus(J)==2)))==0)j=j+1;xchoice(J)=j;if sum(x0(J))==0, errmsg='x0 not correct.'; return; end;end;end;end;end;end;errx=optimget(options2,'TolX');errcon=optimget(options2,'TolCon');fail=0;c=0;% STEP 2 TERMINIATIONwhile stacksize>0c=c+1;% STEP 3 LOADING OF CSPx0=stackx0(:,stacksize);xlb=stackxlb(:,stacksize);xub=stackxub(:,stacksize);x0(find(x0<xlb))=xlb(find(x0<xlb));x0(find(x0>xub))=xub(find(x0>xub));stacksize=stacksize-1;% STEP 4 RELAXATION% PHASE 1con=BNBCON(x0,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin{:});if abs(con)>errcon & settings(1)~=0[x1 dummyfeasflag]=fmincon('0',x0,A,B,Aeq,Beq,xlb,xub,nonlcon,options1,varargin{:});if settings(3) & feasflag==0con=BNBCON(x1,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin{:});if con<errcon, feasflag=1; end;end;else x1=x0; feasflag=1; end;% PHASE 2if feasflag>0[x z convflag]=fmincon(fun,x1,A,B,Aeq,Beq,xlb,xub,nonlcon,options2,varargin{:});if settings(3) & convflag==0con=BNBCON(x,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin{:});if con<errcon, convflag=1; end;end;else convflag=feasflag; end;% STEP 5 FATHOMINGK = find(xstatus==1 & xlb~=xub);separation=1;if convflag<0 | (convflag==0 & settings(2))% FC 1separation=0;elseif z>=z_incumbent & convflag>0% FC 2separation=0;elseif all(abs(round(x(K))-x(K))<errx) & convflag>0% FC 3z_incumbent = z;x_incumbent = x;separation = 0;end;% STEP 6 SELECTIONif separation == 1 & ~isempty(K)dzsep=-1;for i=1:size(K,1)dxsepc = abs(round(x(K(i)))-x(K(i)));if dxsepc>=errx | convflag==0xsepc = x; xsepc(K(i))=round(x(K(i)));dzsepc = abs(feval(fun,xsepc,varargin{:})-z);if dzsepc>dzsepdzsep=dzsepc;ixsep=K(i);end;end;end;% STEP 7 SEPARATIONif xchoice(ixsep)==0% XCHOICE==0branch=1;domain=[xlb(ixsep) xub(ixsep)];while branch==1xboundary=(domain(1)+domain(2))/2;if x(ixsep)<xboundarydomainA=[domain(1) floor(xboundary)];domainB=[floor(xboundary+1) domain(2)];elsedomainA=[floor(xboundary+1) domain(2)];domainB=[domain(1) floor(xboundary)];end;stacksize=stacksize+1;stackx0(:,stacksize)=x;stackxlb(:,stacksize)=xlb;stackxlb(ixsep,stacksize)=domainB(1);stackxub(:,stacksize)=xub;stackxub(ixsep,stacksize)=domainB(2);if domainA(1)==domainA(2)stacksize=stacksize+1;stackx0(:,stacksize)=x;stackxlb(:,stacksize)=xlb;stackxlb(ixsep,stacksize)=domainA(1);stackxub(:,stacksize)=xub;stackxub(ixsep,stacksize)=domainA(2);branch=0;elsedomain=domainA;branch=1;end;end;else% XCHOICE~=0L=find(xchoice==xchoice(ixsep));M=intersect(K,L);[dummy,N]=sort(x(M));part1=M(N(1:floor(size(N)/2))); part2=M(N(floor(size(N)/2)+1:size(N))); stacksize=stacksize+1;stackx0(:,stacksize)=x;O = (1-sum(stackx0(part1,stacksize)))/size(part1,1);stackx0(part1,stacksize)=stackx0(part1,stacksize)+O;stackxlb(:,stacksize)=xlb;stackxub(:,stacksize)=xub;stackxub(part2,stacksize)=0;stacksize=stacksize+1;stackx0(:,stacksize)=x;O = (1-sum(stackx0(part2,stacksize)))/size(part2,1);stackx0(part2,stacksize)=stackx0(part2,stacksize)+O;stackxlb(:,stacksize)=xlb;stackxub(:,stacksize)=xub;stackxub(part1,stacksize)=0;if size(part2,1)==1, stackxlb(part2,stacksize)=1; end;end;elseif separation==1 & isempty(K)fail=fail+1;end;end;% STEP 8 OUTPUTt=toc;Z = z_incumbent;X = x_incumbent;errmsg='';eval(['warning ',currentwarningstate]);function CON=BNBCON(x,A,B,Aeq,Beq,xlb,xub,nonlcon,varargin); if isempty(A), CON1=[]; else CON1 = max(A*x-B,0); end;if isempty(Aeq), CON2=[]; else CON2 = abs(Aeq*x-Beq); end; CON3 = max(xlb-x,0);CON4 = max(x-xub,0);if isempty(nonlcon)CON5=[]; CON6=[];else[C Ceq]=feval(nonlcon,x,varargin{:});CON5 = max(C,0);CON6 = abs(Ceq);end;CON = max([CON1; CON2; CON3; CON4; CON5; CON6]);。