多目标规划MATLAB2012年wgx

MATLAB求解多目标规划江西师范大学数信学院吴根秀一、0-1规划的MATLAB求解数学模型：MINf’xS.T.Ax=bAeqx=beqx=0,1命令格式：x=bintprog(f)x=bintprog(f,A,b)x=bintprog(f,A,b,Aeq,beq)x=bintprog(f,A,b,Aeq,beq,x0)x=bintprog(f,A,b,Aeq,beq,x0,options)[x,fval]=bintprog(...)[x,fval,exitflag]=bintprog(...)[x,fval,exitflag,output]=bintprog(...)数学模型：MINlambdaS.T.F(x)-weight*lambda=goal(达到目标）Ax=b（线性不等式约束）Aeqx=beq（线性等式约束）C(x)=0（非线性不等式约束）Ceq(x)=0（非线性等式约束）lb=x=ubF=[f1(x),f2(x),…]为多目标的目标函数；F与[C(x),Ceq(x)]都是通过function来定义；命令格式：x=fgoalattain(fun,x0,goal,weight)x=fgoalattain(fun,x0,goal,weight,A,b)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)二、多目标规划的MATLAB求解命令格式：x=fgoalattain(fun,x0,goal,weight)x=fgoalattain(fun,x0,goal,weight,A,b)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,lb,ub,nonlcon)x=fgoalattain(fun,x0,goal,weight,A,b,Aeq,beq,...lb,ub,nonlcon,options)[x,fval]=fgoalattain(...)[x,fval,attainfactor]=fgoalattain(...)[x,fval,attainfactor,exitflag]=fgoalattain(...)[x,fval,attainfactor,exitflag,output]=fgoalattain(...)[x,fval,attainfactor,exitflag,output,lambda]=fgoalattain(...)二、多目标规划的MATLAB求解x=fgoalattain(@myfun,x0,goal,weight,A,b,Aeq,beq,...lb,ub,@mycon)wheremyconisaMATLABfunctionsuchasfunction[c,ceq]=mycon(x)c=...%computenonlinearinequalitiesatx.ceq=...%computenonlinearequalitiesat二、多目标规划的MATLAB求解x=fgoalattain(@myfun,x0,goal,weight)wheremyfunisaMATLABfunctionsuchasfunctionF=myfun(x)F=...%Computefunctionvaluesatx.有关优化参数设置：options=optimset(‘GradObj’,‘on’)目标函数的梯度方向参数设置为‘on’时，用下列函数定义：function[F,G]=myfun(x)F=...%Computethefunctionvaluesatxifnargout1%TwooutputargumentsG=...%GradientsevaluatedatxEndThegradientconsistsofthepartialderivativedF/dxofeachFatthepointx.二、多目标规划的MATLAB求解二、多目标规划的MATLAB求解有关优化参数设置：options=optimset(‘GradConstr’,‘on’）约束条件的梯度方向参数设置为‘on’时，用下列函数定义：function[c,ceq,GC,GCeq]=mycon(x)c=...%Nonlinearinequalitiesatxceq=...%Nonlinearequalitiesatxifnargout2%Nonlconcalledwith4outputsGC=...%GradientsoftheinequalitiesGCeq=...%GradientsoftheequalitiesEnd注意：一般weight=abs(goal)模型：x’=(A+BKC)x+Bu，设计K满足目标：Y=Cx1）循环系统的特征值（由命令eig(A+B*K*C)确定）的目标为goal=[-5,-3,-1]2）K中元素均在[-4,4]中;设特征值的weight=abs(goal),定义目标函数F如下：functionF=eigfun(K,A,B,C)F=sort(eig(A+B*K*C));%Evaluateobjectives，由小到大排列优化程序为：A=[-0.500;0-210;01-2];B=[10;-22;01];C=[100;001];K0=[-1-1;-1-1];%Initializecontrollermatrixgoal=[-5-3-1];%Setgoalvaluesfortheeigenvaluesweight=abs(goal)%Setweightforsamepercentagelb=-4*ones(size(K0));%Setlowerboundsonthecontrollerub=4*ones(size(K0));%Setupperboundsonthecontrolleroptions=optimset('Display','iter');%Setdisplayparameter[K,fval,attainfactor]=fgoalattain(@(K)eigfun(K,A,B,C)...goal,weight,[],[],[],[],lb,ub,[],options)二、举例---有关循环控制系统优化问题运行结果如下Activeconstraints:124910K=-4.0000-0.2564-4.0000-4.0000fval=-6.9313-4.1588-1.4099attainfactor=-0.3863二、举例---有关循环控制系统优化问题如果至少保证38.63%的目标精确匹配，设置‘GoalsExactAchieve’参数值为3options=optimset('GoalsExactAchieve',3);[K,fval,attainfactor]=fgoalattain(...@(K)eigfun(K,A,B,C),K0,goal,weight,[],[],[],[],lb,ub,[],...options)Afteraboutseveniterations,asolutionisK=-1.59541.2040-0.4201-2.9046fval=-5.0000-3.0000-1.0000attainfactor=1.0859e-20表明目标已完全匹配二、举例---有关循环控制系统优化问题初等模型举例常见类型定性模型经验公式（拟合、插值）量纲分析比例模型§2.1崖高的估算假如你站在崖顶且身上带着一只具有跑表功能的计算器，你也许会出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度，假定你能准确地测定时间，你又怎样来推算山崖的高度呢，请你分析一下这一问题。我有一只具有跑表功能的计算器。方法一假定空气阻力不计，可以直接利用自由落体运动的公式来计算。例如，设t=4秒，g=9.81米/秒2，则可求得h≈78.5米。221gth我学过微积分，我可以做得更好，呵呵。vKmgdtdvmF除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，此时可设空气阻力正比于石块下落的速度，阻力系数K为常数，因而，由牛顿第二定律可得：kgcevkt令k=K/m，解得代入初始条件v(0)=0，得c=－g/k，故有ktekgkgv再积分一次，得：cekgtkghkt2若设k=0.05并仍设t=4秒，则可求得h≈73.6米。听到回声再按跑表，计算得到的时间中包含了反应时间进一步深入考虑不妨设平均反应时间为0.1秒，假如仍设t=4秒，扣除反应时间后应为3.9秒，代入式①，求得h≈69.9米。222)1(kgektkgkgekgtkghktkt①多测几次，取平均值再一步深入考虑代入初始条件h(0)=0，得到计算山崖高度的公式：将e-kt用泰勒公式展开并令k→0+，即可得出前面不考虑空气阻力时的结果。还应考虑回声传回来所需要的时间。为此，令石块下落的真正时间为t1，声音传回来的时间记为t2，还得解一个方程组：933401212211.ttthkg)ekt(kghkt这一方程组是非线性的，求解不太容易，为了估算崖高竟要去解一个非线性主程组似乎不合情理相对于石块速度，声音速度要快得多，我们可用方法二先求一次h，令t2=h/340，校正t，求石块下落时间t1≈t-t2将t1代入式①再算一次，得出崖高的近似值。例如，若h=69.9米，则t2≈0.21秒，故t1≈3.69秒，求得h≈62.3米。§2.2录像带还能录多长时间录像机上有一个四位计数器，一盘180分钟的录像带在开始计数时为0000，到结束时计数为1849，实际走时为185分20秒。我们从0084观察到0147共用时间3分21秒。若录像机目前的计数为1428，问是否还能录下一个60分钟的节目？rθRl由Wvt)rπ(R22得到212rπvtWR又因和得θRlΔΔtvlΔΔtRvθΔΔ积分得到tθdt)rπWvtv(dθ02120rrπWvtW2πrπWvtW2πθt2212021)()(即从而有rrπWvtWπθn212)(12我们希望建立一个录像带已录像时间t与计数器计数n之间的函数关系。为建立一个正确的模型，首先必须搞清哪些量是常量，哪些量是变量。首先，录像带的厚度W是常量，它被绕在一个半径为r的园盘上，见图。磁带转动中线速度v显然也是常数，否则图象声音必然会失真。此外，计数器的读数n与转过的圈数有关，从而与转过的角度θ成正比。rθRlrrπWvtWπθn212)(12此式中的三个参数W、v和r均不易精确测得，虽然我们可以从上式解出t与n的函数关系，但效果不佳，故令则可将上式简化为：πWvαvW/πrβ2ββtαn故nαφnαββαnt21222令21αaαβb2上式又可化简记成t=an2+bnt=an2+bnrθRl上式以a、b为参数显然是一个十分明智的做法，它为公式的最终确立即参数求解提供了方便。将已知条件代入，得方程组：3.351471478484185.331849(1849)12122tbatbaba从后两式中消去t1，解得a=0.0000291,b=0.04646,故t=0.0000291n2+0.04646n，令n=1428，得到t=125.69（分）由于一盒录像带实际可录像时间为185.33分，故尚可录像时间为59.64分，已不能再录下一个60分钟的