用matlab求解优化问题.doc

8.1.1 线性规划问题的MATLAB求解方法与一般线性规划理论一样，在MATLAB中有线性规划的标准型。在调用MATLAB线性规划函数linprog时，要遵循MATLAB中对标准性的要求。线性规划问题的MATLAB标准形为：在上述模型中，有一个需要极小化的目标函数，以及需要满足的约束条件假设为维设计变量，且线性规划问题具有不等式约束个，等式约束个，那么： 和均为维列向量，为维列向量，为m2维列向量，为维矩阵，为维矩阵需要注意的是：MATLAB标准型是对目标函数求极小，如果遇到是对目标函数求极大的问题，在使用MATLAB求解时，需要在函数前面加一个负号转化为对目标函数求极小的问题；MATLAB标准型中的不等式约束形式为，如果在线性规划问题中出现形式的不等式约束，则我们需要在两边乘以(-1)使其转化为MATLAB中的形式。如果在线性规划问题中出现了“”的约束形式，则我们需要通过添加松弛变量使得不等式约束变为等式约束之后，我们只需要将所有的约束（包括不等式约束和等式约束）转化为矩阵形式的即可。例如，对于如下线性规划模型：要转化为MATLAB标准形，则要经过：（1）原问题是对目标函数求极大，故添加负号使目标变为：；（2）原问题中存在“”的约束条件，故添加负号使其变为： 用MATLAB表达则为c=-4; 2; -1;%将目标函数转化为求极小A=2 -1 1; 8 -2 2; b=12; -8; %不等式约束系数矩阵Aeq=-2 0 1; 1 1 0;beq=3; 7;%等式约束系数矩阵lb=0; 0; 0;ub=Inf; Inf; Inf %对设计变量的边界约束MATLAB优化工具箱中求解线性规划问题的命令为linprog，其函数调用方法有多种形式如下所示：x = linprog(c,A,b)x = linprog(c,A,b,Aeq,beq)x = linprog(c,A,b,Aeq,beq,lb,ub)x = linprog(c,A,b,Aeq,beq,lb,ub,x0)x = linprog(c,A,b,Aeq,beq,lb,ub,x0,options)x = linprog(problem)x,fval = linprog(.)x,fval,exitflag = linprog(.)x,fval,exitflag,output = linprog(.)x,fval,exitflag,output,lambda = linprog(.)输入参数MATLAB工具箱中的linprog函数在求解线性规划问题时，提供的参数为：模型参数、初始解参数和算法控制参数。模型参数x、c、lb、ub、b、beq、A和Aeq在MATLAB标准型中已经介绍了其具体物理意义和获得方法x0为线性规划问题的初始解，该设置仅在中型规模算法中有效，而在默认的大型规模算法和单纯形算法中，MATLAB将忽略一切初始解。options为包含算法控制参数的结构变量，我们可以通过optimset命令对这些具体的控制参数进行设置，例如下述格式 options = optimset(param1,value1,param2,value2,.) 该命令格式创建一组控制参数结构变量options，将参数的具体值赋给单引号之间的参数，任何未被指定的参数将被赋值为，参数值为的具体的含义是将该组控制参数传递给优化函数时将使用MATLAB提供的默认值 在线性规划问题中可以用到的设置参数如下表所示表1 linprog中的控制参数设置参数名称参数设置Diagnostics设置是否显示函数优化中的诊断信息，可以选择on或者off（默认值），该功能主要显示一些退出信息，即linprog函数运算终止的原因Display设置显示信息的级别，当该参数值为off时 ，不显示任何输出信息；当参数值为iter时，将显示每一步迭代的输出信息，iter参数值仅对大型规模算法和中型规模的单纯形算法有效；当参数值为final时，仅显示最终的输出信息Simplex当该参数值为on时，函数采用单纯形算法LargeScale设置是否采用大型规模算法，当参数值为on（默认值）时，使用大型规模算法；当参数值为off时，使用中型规模算法MaxIter算法运行中的最大迭代次数，对于大型规模算法，默认值为85，对于单纯形算法，其默认值为10*设计变量的个数，对于中型有效集算法为10*max(设计变量的个数，不等式约束的个数+边界约束的个数)TolFun函数计算终止的误差限，对于大型规模算法其默认值为1e-8，该控制参数对于中型规模的有效集算法无效。输出参数 linprog函数返回的输出参数有x、fval、exitflag、lambda和output。输出参数x为线性规划问题的最优解输出参数fval为线性规划问题在最优解x处的函数值输出参数exitflag返回的是优化函数计算终止时的状态指示，说明算法终止的原因，其取值和其代表的具体原因如下表所示 表2 参数exitflag的取值和物理意义值物理意义1已经收敛到解x0已经达到最大迭代次数限制options.MaxIter-2没有找到问题的可行点-3问题无有限最优解-4在算法执行过程中遇到了NaN值-5原线性规划问题和其对偶问题均不可行-7搜索方向变化太小，无法进一步获得更优解，说明原线性规划问题或者约束条件是病态的输出参数output是一个返回优化过程中相关信息的结构变量，其属性如下表所示表3 参数output所包含信息属性名称属性含义output.iterations优化过程的实际迭代次数output.algorithm优化过程中所采用的具体算法output.cgiterations0（仅用于大型规模算法，为了后向兼容性而设置的参数）output.message退出信息输出参数lambda是返回线性规划问题最优解x处的拉格朗日乘子的一个结构变量,其总维数等于约束条件的个数，其非零分量对应于起作用的约束条件，其属性如表所示。 表4 参数lambda所包含信息属性名称属性含义ineqlin线性不等式约束的拉格朗日乘子eqlin线性等式约束的拉格朗日乘子upper上界约束的拉格朗日乘子lower下界约束的拉格朗日乘子命令详解1）x = linprog(c,A,b) ；该函数调用格式求解线性规划问题：2）x = linprog(c,A,b,Aeq,beq) ；该函数调用格式求解线性规划问题：即该函数调用格式解决的是既含有线性等式约束，又含有线性不等式约束的线性规划问题，如果在线性规划问题中无线性不等式约束，则可以设A=以及b= 3）x = linprog(c,A,b,Aeq,beq,lb,ub)；该函数调用格式求解线性规划问题： 即在线性规划问题的求解过程中进一步考虑了对设计变量的约束，其中lb和ub均是和设计变量维数相同的列向量，如果对设计变量没有上界约束，可以设置ub(i) = Inf，如果没有下界约束则可以设置lb(i) = -Inf，和(2)类似，如果问题中没有等式约束，则可以设Aeq=以及beq=4）x = linprog(c,A,b,Aeq,beq,lb,ub,x0) 在前面调用方法的基础上设置线性规划问题求解的初始解为x0，该参数仅在使用有效集算法时生效，否则当使用默认的内点算法时，将忽略任何初始点，即参数无效。5）x = linprog(c,A,b,Aeq,beq,lb,ub,x0,options) 用options指定的优化参数进行最小化。可以使用optimset来设置这些参数上面的函数调用格式仅设置了最优解这一输出参数，如果需要更多的输出参数，则可以参照下面的调用格式：1）x,fval = linprog(.) 在优化计算结束之时返回线性规划问题在解x处的目标函数值fval。2）x,fval,exitflag = linprog(.) 在优化计算结束之时返回exitflag值，描述函数计算的退出条件。 3）x,fval,exitflag,output = linprog(.) 在优化计算结束之时返回返回结构变量output4）x,fval,exitflag,output,lambda = linprog(.) 在优化计算结束之时返回线性规划问题最优解x处的拉格朗日乘子lambda例1 用MATLAB求解线性规划问题算法：c=-1;-1;%目标函数，为转化为极小，故取目标函数中设计变量的相反数A=1 -2;1 2;%线性不等式约束b=4;8;lb=0;0;%设计变量的边界约束，由于无上界，故设置ub=Inf;Inf; x,fval=linprog(c,A,b,lb,ub)结果：Optimization terminated.x = 6.0000 1.0000fval = -7.0000例2 用MATLAB求解线性规划问题算法：c=-4;-3; %目标函数，为转化为极小，故取目标函数中设计变量的相反数A=3 4;3 3;4 2;%线性不等式约束b=12;10;8;lb=0;0; %设计变量的边界约束，由于无上界，故设置ub=Inf;Infub=Inf;Inf;x,fval,exitflag=linprog(c,A,b,lb,ub)结果：x = 0.8000 2.4000fval = -10.4000exitflag = 1例3 用MATLAB求解线性规划问题算法：c=-1;-3;1; %目标函数，为转化为极小，故取目标函数中设计变量的相反数Aeq=1 1 2;-1 2 1;%线性等式约束beq=4;4;lb=0;0;0;%设计变量的边界约束，由于无上界，故设置ub=Inf;Inf;Infx,fval,exitflag=linprog(c,Aeq,beq,lb,ub)结果：x = 1.3333 2.6667 0.0000fval = -9.3333exitflag = 1例4 用MATLAB求解线性规划问题算法：c=-3;1;1;%目标函数，为转化为极小，故取目标函数中设计变量的相反数A=1 -2 1;4 -1 -2;%线性不等式约束b=11;-3;Aeq=-2 0 1;%线性等式约束beq=1;lb=0;0;0;%设计变量的边界约束，由于无上界，故设置ub=Inf;Inf;Infub=Inf;Inf;Inf;x,fval,exitflag,output,lambda=linprog(c,A,b,Aeq,beq,lb,ub)结果：x = 4.0000 1.0000 9.0000fval = -2.0000exitflag = 1output = %输出信息结构变量 iterations: 6 %说明优化算法迭代6次 algorithm: large-scale: interior point %说明采用的是大型规模的内点算法 cgiterations: 0 %取值为0，为了向后兼容而设定 message: Optimization terminated. %退出信息lambda = %最优解x处的拉格朗日乘子结构变量 ineqlin: 2x1 double eqlin: -0.6667 upper: 3x1 double lower: 3x1 double8.1.2 整数规划的MATLAB求解方法（一） 用MATLAB求解一般混合整数规划问题由于MATLAB优化工具箱中并未提供求解纯整数规划和混合整数规划的函数，因而需要自行根据需要和设定相关的算法来实现。现在有许多用户发布的工具箱可以解决该类问题，例如比较著名的YALMIP，读者可以自行到网上下载相关的工具包并进行学习。这里我们给出开罗大学的Sherif和Tawfik在MATLAB Central上发布的一个用于求解一般混合整数规划的程序，在此命名为intprog，笔者在原程序的基础上做了简单的修改，将其选择分枝变量的算法由自然序改造成分枝变量选择原则中的一种，即：选择与整数值相差最大的非整数变量首先进行分枝。intprog函数的调用格式如下： x,fval,exitflag=intprog(c,A,b,Aeq,beq,lb,ub,M,TolXInteger) 该函数解决的整数规划问题为：在上述标准问题中，假设为维设计变量，且问题具有不等式约束个，等式约束个，那么：、均为维列向量，为维列向量，为维列向量，为维矩阵，为维矩阵。在该函数中，输入参数有c,A,b,Aeq,beq,lb,ub,M和TolXInteger。其中c为目标函数所对应设计变量的系数，A为不等式约束条件方程组构成的系数矩阵，b为不等式约束条件方程组右边的值构成的向量。Aeq为等式约束方程组构成的系数矩阵，beq为等式约束条件方程组右边的值构成的向量。lb和ub为设计变量对应的上界和下界。M为具有整数约束条件限制的设计变量的序号，例如问题中设计变量为，要求和为整数，则M=2;3;6；若要求全为整数，则M=1:6，或者M=1;2;3;4;5;6。TolXInteger为判定整数的误差限，即若某数x和最邻近整数相差小于该误差限，则认为x即为该整数。在该函数中，输出参数有x, fval和exitflag。其中x为整数规划问题的最优解向量，fval为整数规划问题的目标函数在最优解向量x处的函数值，exitflag为函数计算终止时的状态指示变量。在MATLAB中实现intprog的代码和分析如下：%整数规划的MATLAB实现%Originally Designed By Sherif A. Tawfik,Revised By LiMing, 2009-12-29%函数调用形式x,fval,exitflag=intprog(f,A,b,Aeq,beq,lb,ub,M,TolXInteger)%函数求解如下形式的整数规划问题%min f*x%subject to% A*x=b% Aeq*x=beq% lb=x=ub% M存储有整数约束的变量编号的向量% TolXInteger是判定整数的误差限%函数返回变量%x:整数规划的最优解%fval:目标函数在最优解处的函数值%exitflag =1 收敛到解x% 0 达到线性规划的最大迭代次数% -1 无解%function x,fval,exitflag=intprog(f,A,b,Aeq,beq,lb,ub,M,TolXInteger)%设置不显示求解线性规划过程中的提示信息options = optimset(display,off);%上界的初始值bound=inf; %求解原问题P0的松弛线性规划Q0，首先获得问题的初始解x0,fval0=linprog(f,A,b,Aeq,beq,lb,ub,options); %利用递归法进行二叉树的遍历，实现分枝定界法对整数规划的求解x,fval,exitflag,b=rec_BranchBound(f,A,b,Aeq,beq,lb,ub,x0,fval0,M,TolXInteger,bound);%分枝定界法的递归算法,x为问题的初始解，v是目标函数在x处的取值function xx,fval,exitflag,bb=rec_BranchBound(f,A,b,Aeq,beq,lb,ub,x,v,M,TolXInteger,bound)options = optimset(display,off);%求解不考虑整数约束的松弛线性规划x0,fval0,exitflag0=linprog(f,A,b,Aeq,beq,lb,ub,options);%如果算法结束状态指示变量为负值，即表示无可行解，返回初始输入%或者所目标函数值大于已经获得的上界，返回初始输入if exitflag0bound xx=x; fval=v; exitflag=exitflag0; bb=bound; return;end%确定所有变量是否均为整数，如是，则返回ind=find(abs(x0(M)-round(x0(M)TolXInteger); %该条件表示x0(M)不是整数%如果都是整数if isempty(ind) exitflag=1;%如果当前的解优于已知的最优解，则将当前解作为最优解 if fval0flag br_var=tempbr_var; br_value=tempbr_value; flag=temp_flag; endendif isempty(A) r c=size(Aeq);else r c=size(A);end%分枝后第1个子问题的设置,%添加约束条件Xi=ceil(Xi) ，i即为上面找到的最接近整数的非整数设计变量的序号 A2=A;zeros(1,c);A2(end,br_var)=-1;b2=b;-ceil(br_value); %分枝后的第一个子问题的递归求解x1,fval1,exitflag1,bound1=rec_BranchBound(f,A1,b1,Aeq,beq,lb,ub,x0,fval0,M,TolXInteger,bound);exitflag=exitflag1;if exitflag10 & bound10 & bound2 distance(hopS) + transmat(hopS,hopE) distance(hopE) = distance(hopS) + transmat(hopS,hopE); parent(hopE) = hopS; queue = queue hopE; end endend% 回溯进行最短路径的查找r_path = pathE; i = parent(pathE);while i=pathS & i=0 r_path = i r_path; i = parent(i);end if i=pathS r_path = i r_path;else r_path = end% 返回最短路径的权和r_cost = distance(pathE);例1旅行家问题假设有一个旅行家，他居住在城市，他计划游览他附近的若干城市，假设这些城市网络可以用下图来表示，现在他首先要估算从A出发到各城市的最少花费，然后确定行程，我们如何利用Dijkstra算法来解决他的问题呢？即如何求解到各顶点的最短距离。解：算法：W=0 1 2 Inf 7 4 8 Inf ; 1 0 2 5 Inf Inf 7 Inf ; 2 2 0 1 5 Inf Inf Inf; Inf 5 1 0 3 Inf Inf 6 ; 7 Inf 5 3 0 3 Inf 4 ; 4 Inf Inf Inf 3 0 2 6; 8 7 Inf Inf Inf 2 0 4; Inf Inf Inf 6 4 6 4 0;r_path, r_cost = dijkstra(1, 8, W)结果：r_path = 1 3 4 8r_cost = 9与理论分析结果相同。2、Floyd算法根据前面的算法描述，可以用MATLAB编写程序如下，M-函数文件floydSPR.m的输入参数为图的邻接矩阵，输出参数即为算法描述中的距离矩阵D和最短路径矩阵R。%floydSPR算法%a-赋权邻接矩阵%D-为距离矩阵%R-最短路径矩阵%By GreenSim Groupfunction D,R=floydSPR(a)n=size(a,1);D=a;R=zeros(n,n);for i=1:n for j=1:n if D(i,j)=inf R(i,j)=j; end endendfor k=1:n for i=1:n for j=1:n if D(i,k)+D(k,j)0&f(i,j)=0) a(i,j)=b(i,j); elseif(C(i,j)0&f(i,j)=C(i,j) a(j