您的位置：360文档中心› 大连理工优化方法大作业MATLAB编程

大连理工优化方法大作业MATLAB编程

合集下载

最优化方法及其matlab程序设计

最优化方法及其matlab程序设计

最优化方法及其matlab程序设计
最优化方法是一种利用各种技术，以提高某项工作，工程或系统
的效率为目标，并让其在某些给定基准测试中改善性能的过程。

它可
以用来提高计算机系统的性能，减少加工时间，提高生产率，等等。

Matlab是一种非常适用于最优化的程序设计语言，它拥有许多强
大的分析功能，例如数值分析、线性规划、非线性规划、二次规划、
优化算法、深度学习、图形处理和仿真等。

因此，Matlab可以帮助用
户找到最优解决方案，比如解决所谓的NP难问题，这些问题很难在
“合理”时间内找到最优解。

要在matlab中实现最优化方法，首先要定义和描述优化问题。

然后，选择合适的优化器。

一般来说，FMINCON函数可以满足大多数最优
化问题的要求，因为它可以通过求解约束和非线性问题来实现最优化。

在函数中，用户可以指定具体的约束条件、目标函数、初始解和其他
一些参数，以便更好地进行最优化。

此外，matlab中还提供了其他一些有用的优化函数，可以用于解
决更复杂的问题，包括FMINUNC、FMINBND等。

这些函数都可以实现更
高级的最优化算法，例如迭代算法、模拟退火算法、遗传算法等。

最后，用户还可以使用matlab自带的toolbox来进行最优化，例
如Optimization Toolbox。

这个工具包可以帮助用户调整参数，从而
实现最优解。

同时，它还提供了有关具体优化策略的解释，以便了解
该策略的实现方法以及它的应用范围。

总的来说，matlab可以实现各种最优化方法，无论是简单的还是
复杂的，都可以通过它找到最佳解决方案。

2016年大连理工大学优化方法上机大作业

2016年大连理工大学优化方法上机大作业

2016年理工大学优化方法上机大作业学院：专业：班级：学号：：上机大作业1：1.最速下降法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = steepest(x0,eps)gk = grad(x0);res = norm(gk);k = 0;while res > eps && k<=1000dk = -gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;gk = grad(xk);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x=steepest(x0,eps)2.牛顿法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad2(x)g = zeros(2,2);g(1,1)=2+400*(3*x(1)^2-x(2));g(1,2)=-400*x(1);g(2,1)=-400*x(1);g(2,2)=200;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = newton(x0,eps)gk = grad(x0);bk = [grad2(x0)]^(-1);res = norm(gk);k = 0;while res > eps && k<=1000dk=-bk*gk;xk=x0+dk;k = k+1;x0 = xk;gk = grad(xk);bk = [grad2(xk)]^(-1);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x1=newton(x0,eps)--The 1-th iter, the residual is 447.213595--The 2-th iter, the residual is 0.000000x1 =1.00001.00003.BFGS法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = bfgs(x0,eps)g0 = grad(x0);gk=g0;res = norm(gk);Hk=eye(2);k = 0;while res > eps && k<=1000dk = -Hk*gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;fa0=xk-x0;x0 = xk;go=gk;gk = grad(xk);y0=gk-g0;Hk=((eye(2)-fa0*(y0)')/((fa0)'*(y0)))*((eye(2)-(y0)*(fa0)')/((fa0)'*( y0)))+(fa0*(fa0)')/((fa0)'*(y0));res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;End>> clear>> x0=[0,0]';>> eps=1e-4;>> x=bfgs(x0,eps)4.共轭梯度法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star =CG(x0,eps)gk = grad(x0);res = norm(gk);k = 0;dk = -gk;while res > eps && k<=1000ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;g0=gk;gk = grad(xk);res = norm(gk);p=(gk/g0)^2;dk1=dk;dk=-gk+p*dk1;fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk; end>> clear>> x0=[0,0]'; >> eps=1e-4; >> x=CG(x0,eps)上机大作业2：function f= obj(x)f=4*x(1)-x(2)^2-12;endfunction [h,g] =constrains(x) h=x(1)^2+x(2)^2-25;g=zeros(3,1);g(1)=-10*x(1)+x(1)^2-10*x(2)+x(2)^2+34;g(2)=-x(1);g(3)=-x(2);endfunction f=alobj(x) %拉格朗日增广函数%N_equ等式约束个数?%N_inequ不等式约束个数N_equ=1;N_inequ=3;global r_al pena;%全局变量h_equ=0;h_inequ=0;[h,g]=constrains(x);%等式约束部分?for i=1:N_equh_equ=h_equ+h(i)*r_al(i)+(pena/2)*h(i).^2;end%不等式约束部分for i=1:N_inequh_inequ=h_inequ+(0.5/pena)*(max(0,(r_al(i)+pena*g(i))).^2-r_al(i).^2) ;end%拉格朗日增广函数值f=obj(x)+h_equ+h_inequ;function f=compare(x)global r_al pena N_equ N_inequ;N_equ=1;N_inequ=3;h_inequ=zeros(3,1);[h,g]=constrains(x);%等式部分for i=1:1h_equ=abs(h(i));end%不等式部分for i=1:3h_inequ=abs(max(g(i),-r_al(i+1)/pena));endh1 = max(h_inequ);f= max(abs(h_equ),h1); %sqrt(h_equ+h_inequ);function [ x,fmin,k] =almain(x_al)%本程序为拉格朗日乘子算法示例算法%函数输入：% x_al:初始迭代点% r_al：初始拉格朗日乘子N-equ:等式约束个数N_inequ:不等式约束个数?%函数输出% X：最优函数点FVAL:最优函数值%============================程序开始================================ global r_al pena ; %参数（全局变量）pena=10; %惩罚系数r_al=[1,1,1,1];c_scale=2; %乘法系数乘数cta=0.5; %下降标准系数e_al=1e-4; %误差控制围max_itera=25;out_itera=1; %迭代次数%===========================算法迭代开始============================= while out_itera<max_iterax_al0=x_al;r_al0=r_al;%判断函数?compareFlag=compare(x_al0);%无约束的拟牛顿法BFGS[X,fmin]=fminunc(alobj,x_al0);x_al=X; %得到新迭代点%判断停止条件?if compare(x_al)<e_aldisp('we get the opt point');breakend%c判断函数下降度?if compare(x_al)<cta*compareFlagpena=1*pena; %可以根据需要修改惩罚系数变量elsepena=min(1000,c_scale*pena); %%乘法系数最大1000disp('pena=2*pena');end%%?更新拉格朗日乘子[h,g]=constrains(x_al);for i=1:1%%等式约束部分r_al(i)= r_al0(i)+pena*h(i);endfor i=1:3%%不等式约束部分r_al(i+1)=max(0,(r_al0(i+1)+pena*g(i)));endout_itera=out_itera+1;end%+++++++++++++++++++++++++++迭代结束+++++++++++++++++++++++++++++++++ disp('the iteration number');k=out_itera;disp('the value of constrains'); compare(x_al)disp('the opt point');x=x_al;fmin=obj(X);>> clear>> x_al=[0,0];>> [x,fmin,k]=almain(x_al)上机大作业3：1、>> clear alln=3; c=[-3,-1,-3]'; A=[2,1,1;1,2,3;2,2,1;-1,0,0;0,-1,0;0,0,-1];b=[2,5,6,0,0,0]'; cvx_beginvariable x(n)minimize( c'*x)subject toA*x<=bcvx_endCalling SDPT3 4.0: 6 variables, 3 equality constraints------------------------------------------------------------num. of constraints = 3dim. of linear var = 6*******************************************************************SDPT3: Infeasible path-following algorithms*******************************************************************version predcorr gam expon scale_dataNT 1 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime-------------------------------------------------------------------0|0.000|0.000|1.1e+01|5.1e+00|6.0e+02|-7.000000e+01 0.000000e+00| 0:0:00| chol 1 11|0.912|1.000|9.4e-01|4.6e-02|6.5e+01|-5.606627e+00 -2.967567e+01| 0:0:01| chol 1 12|1.000|1.000|1.3e-07|4.6e-03|8.5e+00|-2.723981e+00 -1.113509e+01| 0:0:01| chol 1 13|1.000|0.961|2.3e-08|6.2e-04|1.8e+00|-4.348354e+00 -6.122853e+00| 0:0:01| chol 1 14|0.881|1.000|2.2e-08|4.6e-05|3.7e-01|-5.255152e+00 -5.622375e+00| 0:0:01| chol 1 15|0.995|0.962|1.6e-09|6.2e-06|1.5e-02|-5.394782e+00 -5.409213e+00| 0:0:01| chol 1 16|0.989|0.989|2.7e-10|5.2e-07|1.7e-04|-5.399940e+00 -5.400100e+00| 0:0:01| chol 1 17|0.989|0.989|5.3e-11|5.8e-09|1.8e-06|-5.399999e+00 -5.400001e+00| 0:0:01| chol 1 18|1.000|0.994|2.8e-13|4.3e-11|2.7e-08|-5.400000e+00 -5.400000e+00| 0:0:01| stop: max(relative gap, infeasibilities) < 1.49e-08-------------------------------------------------------------------number of iterations = 8primal objective value = -5.39999999e+00dual objective value = -5.40000002e+00gap := trace(XZ) = 2.66e-08relative gap = 2.26e-09actual relative gap = 2.21e-09rel. primal infeas (scaled problem) = 2.77e-13rel. dual " " " = 4.31e-11rel. primal infeas (unscaled problem) = 0.00e+00rel. dual " " " = 0.00e+00norm(X), norm(y), norm(Z) = 4.3e+00, 1.3e+00, 1.9e+00norm(A), norm(b), norm(C) = 6.7e+00, 9.1e+00, 5.4e+00Total CPU time (secs) = 0.71CPU time per iteration = 0.09termination code = 0DIMACS: 3.6e-13 0.0e+00 5.8e-11 0.0e+00 2.2e-09 2.3e-09-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval): -5.42、>> clear alln=2; c=[-2,-4]'; G=[0.5,0;0,1]; A=[1,1;-1,0;0,-1]; b=[1,0,0]'; cvx_beginvariable x(n)minimize( x'*G*x+c'*x)subject toA*x<=bcvx_endCalling SDPT3 4.0: 7 variables, 3 equality constraintsFor improved efficiency, SDPT3 is solving the dual problem.------------------------------------------------------------num. of constraints = 3dim. of socp var = 4, num. of socp blk = 1dim. of linear var = 3******************************************************************* SDPT3: Infeasible path-following algorithms*******************************************************************version predcorr gam expon scale_dataNT 1 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime-------------------------------------------------------------------0|0.000|0.000|8.0e-01|6.5e+00|3.1e+02| 1.000000e+01 0.000000e+00| 0:0:00| chol 1 1 1|1.000|0.987|4.3e-07|1.5e-01|1.6e+01| 9.043148e+00 -2.714056e-01| 0:0:00| chol 1 1 2|1.000|1.000|2.6e-07|7.6e-03|1.4e+00| 1.234938e+00 -5.011630e-02| 0:0:00| chol 1 1 3|1.000|1.000|2.4e-07|7.6e-04|3.0e-01| 4.166959e-01 1.181563e-01| 0:0:00| chol 1 1 4|0.892|0.877|6.4e-08|1.6e-04|5.2e-02| 2.773022e-01 2.265122e-01| 0:0:00| chol 1 1 5|1.000|1.000|1.0e-08|7.6e-06|1.5e-02| 2.579468e-01 2.427203e-01| 0:0:00| chol 1 1 6|0.905|0.904|3.1e-09|1.4e-06|2.3e-03| 2.511936e-01 2.488619e-01| 0:0:00| chol 1 1 7|1.000|1.000|6.1e-09|7.7e-08|6.6e-04| 2.503336e-01 2.496718e-01| 0:0:00| chol 1 1 8|0.903|0.903|1.8e-09|1.5e-08|1.0e-04| 2.500507e-01 2.499497e-01| 0:0:00| chol 1 19|1.000|1.000|4.9e-10|3.5e-10|2.9e-05| 2.500143e-01 2.499857e-01| 0:0:00| chol 1 1 10|0.904|0.904|4.7e-11|1.3e-10|4.4e-06| 2.500022e-01 2.499978e-01| 0:0:00| chol 2 2 11|1.000|1.000|2.3e-12|9.4e-12|1.2e-06| 2.500006e-01 2.499994e-01| 0:0:00| chol 2 2 12|1.000|1.000|4.7e-13|1.0e-12|1.8e-07| 2.500001e-01 2.499999e-01| 0:0:00| chol 2 2 13|1.000|1.000|2.0e-12|1.0e-12|4.2e-08| 2.500000e-01 2.500000e-01| 0:0:00| chol 2 2 14|1.000|1.000|2.6e-12|1.0e-12|7.3e-09| 2.500000e-01 2.500000e-01| 0:0:00|stop: max(relative gap, infeasibilities) < 1.49e-08-------------------------------------------------------------------number of iterations = 14primal objective value = 2.50000004e-01dual objective value = 2.49999996e-01gap := trace(XZ) = 7.29e-09relative gap = 4.86e-09actual relative gap = 4.86e-09rel. primal infeas (scaled problem) = 2.63e-12rel. dual " " " = 1.00e-12rel. primal infeas (unscaled problem) = 0.00e+00rel. dual " " " = 0.00e+00norm(X), norm(y), norm(Z) = 3.2e+00, 1.5e+00, 1.9e+00norm(A), norm(b), norm(C) = 3.9e+00, 4.2e+00, 2.6e+00Total CPU time (secs) = 0.36CPU time per iteration = 0.03termination code = 0DIMACS: 3.7e-12 0.0e+00 1.3e-12 0.0e+00 4.9e-09 4.9e-09------------------------------------------------------------------------------------------------------------------------------- Status: SolvedOptimal value (cvx_optval): -3。

MatLab与Fortran混合编程实现结构优化和可靠性分析

MatLab与Fortran混合编程实现结构优化和可靠性分析
ＺＨＡＯ ｕ ｌ，ＹＩＰ ｎ Ｘｉ—ｉ ｉｇ
（ ａｕｙｏ ｏｓｕｔｎＥ ｇｎｅｎ Ｄ ｌ ｎＵ ｉｒｔ ｏ Ｔｃｎｌ ｙ ａ ａ ，Ｌ ｏｉ １０４ ｈｎ ） Ｆｃｌ ｔ ｆＣｎｔ ｃｏ ｎｉ ｒ ｇ， ａｉ ｎ ｅ ｉ ｅｈｏ ｇ ，Ｄ ｌ ｎ ｉ ｎ ｇ １６２ ，Ｃ ｉ ｒ ｉ ｅｉ ａ ｖ ｓｙｆ ｏ ｉ ａ ｎ ａ
２００ 第 第 ３ １卷 ６月 １２年 翅
Ｊｕ ｌ ｆ ｔ ｅｏｒｅ ａｄＡ ｃｉｃ 珀 ｏｍａ ０ ｗａ水利与建筑工程学报 ｅ Ｒ ｓＩｃｓ ｎ ｒ ｔ 】 ｒ ｌ ｈｅｕ
Ｖ １１ ｏ３ ｏ．０Ｎ ．
Ｊ ｎ ．， ２ ｕ ２０ １
Ｍａ ａ ｔ ｂ与 Ｆ ｒ ａ Ｌ ｏｔ ｎ混 合 编 程 实 现 ｒ 结 构 优 化 和 可 靠 性 分 析
ｓｒ ｔ ｕｃｔ ｅ ｕｒ ｓ．
Ｋ ｅ ｗｏ ｄ ｙ ｒ ｓ：Ｍ ａ Ｌａ ａ ｕ ｇ ｔ ｂ ｌｎｇ ａ ｅ；Ｆｏ ｔａ ａ ｇ ａ ｅ；ｍ ｉｉ ｏ ｒ ｍ ｍ ｉ ｇ； ｓｒ ｃ ｕ ａ ｐ ｉ ｉａ ｉ ｎ；ｓｒ ｔ ａ ｅｉ ｒｒ ｎ ｌ ｎ ｕｇ ｘ ｎｇ ｐｒ ｇ ａ ｎ ｔ ｕ ｔ ｒ ｌｏ ｔｍ ｚ ｔｏ ｔｕｃｕｒ ｌｒ ｌ－
ａ ｉｉ ｎａｙ ｉ ｂ ｌ ｙ ａ ｌ ｓｓ ｔ
Ｍ ｔａ 化 工 具 箱 针 对 各 种 优 化 问题 给 出 了 ａＬｂ优 完整 的解 决方 案 ｌ＿ ， 内容 涵盖 了包括最 值 问题 、 １ 其 ｌ
Ｆｒａ 有 限元 程序成 功 连 接 ； ｏｒ ｔｎ 然后 在 Ｍａ ａ ｔ ｂ环 境 下 ｌ
赵 秀 丽 ， 平 易
（ 大连理工大学 建设工程学部 ， 辽宁 大连 １６２ ） １０４
摘 要 ： 为有效利用 Ｍａ ａ ｔ ｂ和现有 Ｆ ￣ａ 有 限元程 序 ， Ｌ ｏｒ ｎ 采用 Ｍａ ａ 与 Ｆｒａ 合编程 。通过 Ｍａ ａ 程 ｔｂ Ｌ ｏｒ ｔ ｔｎ 昆 ｔｂ Ｌ 序接 口， Ｍａ ａ 与 Ｆ￣ａ 有 限元 程序成功连接实现平 面桁架结 构 的静 力分析 ； 将 ｔｂ Ｌ ｏｒ ｎ 然后 在 Ｍ ｔａ ａ￣ ｌｂ环境下 运用优 化工 具箱 和统计工具箱编程实现 了平面桁架 的结 构优 化设 计和可靠指标的计算。 关键词 ： ａ ａ 语言 ；ｏ ｒ 语 言 ； Ｍｔ ｂ Ｌ Ｆ￣ａ ｎ 混合编程 ； 结构优化 ； 构可靠性 分析 结

(整理)大连理工大学--优化作业----程序.

(整理)大连理工大学--优化作业----程序.

1.1程序（Java）public class Wolfe_Powell {public static double getFx ( double[] x ) {double x1= x[0]; double x2 = x[1];double Fx= 100 * (x2-x1*x1)* (x2-x1*x1) + (1-x1)* (1-x1) ;return Fx;}public static double[] getDeltFx ( double[] x ) {double x1= x[0]; double x2 = x[1];double [] deltFx = new double[2];deltFx [0] = -400*(x2 - x1* x1) *x1- 2*(1- x1) ;deltFx [1] = 200*(x2- x1 * x1) ;return deltFx ;}public static double getDeltFx_Sk ( double[] deltFx , double[] Sk ) { double a = 0 ;for ( int i = 0 ; i < Sk.length ; i++ ) {a = a + deltFx[ i ] * Sk[ i ] ;}return a ;}public static double getL ( double[] x, double[] s ) {double x1= x[0]; double x2 = x[1];double c1 =0.1 , c2 =0.5 ,a =0 , b=1e8 ,L= 1;double Fx0 , Fx1 ，deltFx1_Sk ,deltFx0_Sk ,temp ,temp2;double[] deltFx0 , deltFx1 ;Fx0 = getFx(x) ;deltFx0 = getDeltFx (x) ;deltFx0_Sk = getDeltFx_Sk( deltFx0 , s) ;temp = c2 * getDeltFx_Sk( deltFx0 , s) ;for ( int i=0;i< 1e8 ; i++){temp2 = -c1 * L * deltFx0_Sk ;x[0] = x1 + L *s[0] ;x[1] = x2 + L *s[1] ;Fx1 = getFx(x) ;deltFx1 = getDeltFx (x) ;deltFx1_Sk = getDeltFx_Sk (deltFx1 , s) ;if( (Fx0 - Fx1 ) >= temp2 && deltFx1_Sk >= temp){break ;}else if( (Fx0 - Fx1 ) < temp2 ){b = L ;L = (L +a) /2 ;}else if ( deltFx1_Sk < temp ) {a = L ;L = ( L + b ) / 2 >= 2*L ? (2*L):( L + b ) / 2;}}System.out.println(" L= " + L);System.out.println(" 计算次数" + i );return L ;}public static void main(String[] args) {Wolfe_Powell temp =new Wolfe_Powell();double [] X = { -1 ,1 } ; double [] sk = { 1 ,1} ; temp.getL( X ,sk) ;}}1.2实验结果步长L = 0.00390625 x =[-0.9992 , 1.0324] 计算次数82.1程序（Java）public class GongE {public static double getFx ( double[] x ) {double x1= x[0];double x2 = x[1];double Fx= x1*x1 - 2*x1*x2 + 2*x2*x2 +x3*x3 - x2*x3 +2 * x1 +3*x2 -x3 ; return Fx;}public static double[] getDeltFx ( double[] x ) {double x1= x[0];double x2 = x[1];double [] deltFx = new double[x.length];deltFx [0] = 2*x1 - 2*x2+2 ;deltFx [1] = -2*x1 +4*x2 - x3 +3;deltFx [2] = 2*x3 -x2 -1 ;return deltFx ;}public static double[] getX ( double[] x ) {double[] g0,g1;double[] s0= new double[x.length];double[] s1=new double[x.length];double g0_L,g1_L ,L ,temp;double[] x0 =x ;int k =0 ;g0 = getDeltFx ( x0 ) ;for ( int j = 0 ; j < x.length ; j++ ) {s0[ j ] = -g0[ j ] ;}for (int i = 0 ;i<2; i ++,k++){g0 = getDeltFx ( x0 ) ;g0_L = getDeltFx_Sk ( s0 , s0 ) ;L =getL(x0,s0); // 例题一中的方法取得步长Lfor(int j=0;j<x.length ; j++){x0[j]= x0[j]+ s0[j]*L ;}g1 = getDeltFx(x0) ;g1_L = getDeltFx_Sk ( g1 , g1 );if ( Math.sqrt( g1_L )<= 1e-2 ) {break ;}else{temp = g1_L/ g0_L ;for(int j=0;j<x.length ; j++){s0[j] = -g1[j] + temp * s0[j];}} }return x0;}public static void main(String[] args) {GongE temp =new GongE();double [] x = { 1,1 } ;double[] result = temp.getX(x) ;for ( int i = 0 ; i < x.length ; i++ ) {System.out.println ( "result[" + i + "]=" + result[ i ] ) ;} } }2.2实验结果最优点x*=[-4,-3,-1] 最优解f*=-83.1公用程序（Java）public static double getFx ( double[] x ) { //取得Fx 值double x1= x[0]; double x2 = x[1];double Fx = x1 + 2 * x2 * x2 + Math.exp ( x1 * x1 + x2 * x2 ) ;return Fx ;}public static double[] getDeltFx ( double[] x ) { //取得Fx 的梯度值double x1= x[0]; double x2 = x[1];double[] deltFx = new double[ 2 ] ;deltFx[ 0 ] = 1 + 2 * x1 * Math.exp ( x1 * x1 + x2 * x2 ) ;deltFx[ 1 ] = 4 * x2 + 2 * x2 * Math.exp ( x1 * x1 + x2 * x2 ) ;return deltFx ;}3.2.1最速下降法程序（Java）public class FastWay {public static double[] getX ( double[] x ) {double [] deltF0 = new double[2]; double L =0;for ( int i = 0 ; i < 1e1 ; i++ ) {deltF0 = getDeltFx(x);for(int j=0 ;j <deltF0.length ;j++){ //取得负梯度deltF0[j] = - deltF0[j];}L = getL ( x , deltF0 ) ; // 调用习题1的不精确搜索取得步长Lif ( Math.sqrt ( getDeltFx_Sk ( deltF0 , deltF0 ) ) <= 1e-3 ) {System.out.println ( "最终计算次数" + i ) ;System.out.println("x1=" + x[0]+" x2=" + x[1]);break ;}x[0] = x[0]+ L * deltF0[ 0 ] ; x[1]= x[1]+ L * deltF0[ 1 ] ;}return x;}public static void main ( String[] args ) {FastWay temp = new FastWay () ;double[] x0 = { 2 , 2} ; temp.getX(x0) ;}3.2.2最速下降法结果最优点X*=[-0.4194 0] 最优解f*=0.7729 计算次数count=10 3.3.1牛顿法程序（Java）public static double[] getDeltFx ( double[] x ) {double x1 = x[ 0 ] ; double x2 = x[ 1 ] ;double[] one = new double[ 2 ] ;double exp =Math.exp( Math.pow(x1,2)+Math.pow(x2,2)) ;one[ 0 ] = 1+ 2*x1*exp ; one[ 1 ] = 4* x2 +2*x2*exp ;double[][] two = new double[2][2] ;two[0][0] = 2*exp *(1+2*Math.pow(x1,2)) ;two[1][1] = 2*exp *(1+2*Math.pow(x2,2)) +4 ;double[] deltFx = new double[ 2 ] ;for (int i = 0 ; i < 2 ; i++ ) {deltFx[0] = one[ 0 ]/two[0][0] ;deltFx[1] = one[ 1 ]/two[1][1] ;}return deltFx;}public static void main ( String[] args ) {double[] x = { 1 , 0} ;double[] DeltFx = new double [2] ;for(int i =0 ;i <1e3;i++){DeltFx = getDeltFx(x);x[0] = x[0]- DeltFx[0];x[1] = x[1]- DeltFx[1];if( Math.sqrt( getDeltFx_Sk(DeltFx,DeltFx ) ) <= 1e-4){System.out.println("计算次数为" + i);break ;}}System.out.println(" x1= " +x[0] +" x2= " + x[1] +"\n") ;System.out.println(" Fx= " +getFx(x)) ;}3.3.2牛顿法结果最优点X*=[ -0.4194 , 0] 最优解f*= 0.7729 计算次数count=5 3.4.1 BFGS法程序（matlab）function [x,val,k] = bfgs(fun,gfun,x0)maxk=1000; sigma=0.4; rho=0.55 ; epsion=1e-5;k=0 ; n =length(x0); Bk=eye(n); %Bk=feval('Hess',x0);while (k<maxk)gk=feval(gfun,x0);if(norm(gk)<epsion),break;end;dk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk)oldf=feval(fun,x0)if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x)-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);end;k=k+1; x0=x;endval=feval(fun,x0);3.4.2 BFGS法结果最优点X*=[-0.4194 0] 最优解f*=0.7729 计算次数count=44.1 有效集法（matlab）4.1.1 主程序function[x , Lagrange , exitflag , output]= TwoProg (H,c,Ae,be,Ai,bi,x0)n=length(x0); x=x0; ni=length(bi); ne=length(be); Lagrange =zeros(ne+ni,1); index=ones(ni,1); for(i=1:ni)if(Ai(i,:)*x>bi(i)+1e-9),index(i)=0;endend%算法主程序k=0;while(k<=1e4)%求解子问题Temp=[];if(ne>0),Temp=Ae ; endfor(j=1:ni)if(index(j)>0),Temp=[Temp;Ai(j,:)];endendgk=H*x+c;[m1,n1]=size(Temp);[dk,Lagrange ]=SubPro (H,gk , Temp,zeros(m1,1));if(norm(dk)<= 1.0e-6)y=0.0;if(length(Lagrange )>ne)[y,jk]=min(Lagrange (ne+1:length(Lagrange )));endif(y>=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk)index(i)=0;break;endendendk=k+1;elseexitflag=1;%求步长alpha=1.0;tm=1.0;for(i=1:ni)if((index(i)==0)&(Ai(i,:)*dk<0))tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if(tm1<tm)tm=tm1;ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1),index(ti)=1;endendif(exitflag==0),break;endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;4.1.2 目标函数function f=fun(x)x1=x(1); x2=x(2); f=eval ('x1+2*x2^2+exp(x1^2+x2^2)');4.1.3 子问题函数function[x, Lagrange ]= SubPro (H ,c, Ae, be)[m,n]=size(Ae);ginvH=pinv(H);if(m>0)rb=Ae*ginvH*c+be;Lagrange =pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*Lagrange -c);elsex=-ginvH*c;Lagrange =0;end4.1.4 运行函数H=[2 -2;-2 4];c=[-2 -6]';Ae=[ ];be=[ ];Ai=[1 -2;-0.5 -0.5;1 0;0 1];bi=[-2 -1 0 0]';x0=[0 1 ]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)4.2 有效集法结果内部点初始点x0=[0 0] 最优点X*=[0.8 1.2] 最优解f*=-7.2 迭代次数=10 边界点初始点x0=[1 1] 最优点X*=[0.8 1.2] 最优解f*=-7.2 迭代次数=2 检验点初始点x0=[0 1] 最优点X*=[0.8 1.2] 最优解f*=-7.2 迭代次数=75.1 乘子法程序（matlab）5.1.1 chengZi程序---乘子法主程序function[x,mu,Lagrange ,output]=chengZi(fun,hf,gf,dfun,dhf,dgf,x0)sigma=2.0;count=0;innerCount=0;eta=2.0;θ=0.8;%PHR算法中的实参数θx=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(gi);%选取乘子向量的初始值mu=0.1*ones(l,1);Lagrange =0.1*ones(m,1);btak=10;btaold=10;%用来检验终止条件的两个值while(btak>1e-6&count<1e3 )%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs('Lagr','LagrTiDu',x0,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange ,sigma);innerCount=innerCount+ik;he=feval(hf,x);gi=feval(gf,x);btak=0.0;for(i=1:l),btak=btak+he(i)^2; endfor(i=1:m)temp=min(gi(i),Lagrange (i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);if btak>1e-6if(count>=2&btak>θ*btaold)sigma=eta*sigma;end%更新乘子向量for(i=1:l),mu(i)=mu(i)-sigma*he(i);endfor(i=1:m)Lagrange (i)=max(0.0,Lagrange (i)-sigma*gi(i));endendcount=count+1;btaold=btak;x0=x;endf=feval(fun,x)output.inner_iter=innerCount;output.iter=count;output.bta=btak;output.fval=f;5.1.2 f1程序---目标函数function f=f1(x)f=4*x(1)-x(2)^2-12;5.1.3 h1程序---等式约束function he=h1(x)he=25-x(1)^2-x(2)^2;5.1.4 g1程序---不等式约束function gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;5.1.5 df1程序---目标函数的梯度文件function g=df1(x)g=[4 ,-2.0*x(2)]';5.1.6 dhe程序---等式约束（向量）函数的Jacobi矩阵（转置）function dhe=dh1(x)dhe=[-2*x(1),-2.0*x(2)]';5.1.7 dgi程序---不等式约束（向量）函数的Jacobi矩阵（转置）function dgi=dg1(x)dgi=[10-2*x(1),10-2*x(2);0,1;1,0]';5.1.8 LagrTiDu程序---增广拉格朗日函数的梯度程序function result=LagrTiDu(x,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange ,sigma) result=feval(dfun,x);he=feval(hf,x);gi=feval(gf,x);dhe=feval(dhf,x);dgi=feval(dgf,x);l=length(he);m=length(gi);for(i=1:l)result=result+(sigma*he(i)-mu(i))*dhe(:,i);精品文档endfor(i=1:m)result=result+(sigma*gi(i)-Lagrange (i))*dgi(:,i);end5.1.9 Lagr程序---增广拉格朗日函数程序function result=Lagr(x,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange ,sigma)f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);result=f;s1=0.0;for(i=1:l)result=result-he(i)*mu(i);s1=s1+he(i)^2;endresult=result+0.5*sigma*s1;s2=0.0;for(i=1:m)s3=max(0.0,Lagrange (i)-sigma*gi(i));s2=s2+s3^2-Lagrange (i)^2;endresult=result+s2/(2.0*sigma);5.2 乘子法结果初始点x0=[0 , 0] 最优点X*=[1.0013，4.8987] 最优解f*= -31.9923 等式乘子向量L hu=1.0156 不等式乘子向量Lg=0.75445精品文档。

最优化方法的Matlab实现

最优化方法的Matlab实现

最优化方法的Matlab实现Matlab中使用最优化方法可以使用优化工具箱。

在优化工具箱中，有多种最优化算法可供选择，包括线性规划、非线性规划、约束优化等。

下面将详细介绍如何在Matlab中实现最优化方法。

首先，需要建立一个目标函数。

目标函数是最优化问题的核心，它描述了要优化的变量之间的关系。

例如，我们可以定义一个简单的目标函数：```matlabfunction f = objFun(x)f=(x-2)^2+3;end```以上代码定义了一个目标函数`objFun`，它使用了一个变量`x`，并返回了`f`的值。

在这个例子中，目标函数是`(x-2)^2 + 3`。

接下来，需要选择一个最优化算法。

在Matlab中，有多种最优化算法可供选择，如黄金分割法、割线法、牛顿法等。

以下是一个使用黄金分割法的示例：```matlabx0=0;%初始点options = optimset('fminsearch'); % 设定优化选项```除了黄金分割法，还有其他最优化算法可供选择。

例如，可以使用`fminunc`函数调用一个无约束优化算法，或者使用`fmincon`函数调用带约束的优化算法。

对于非线性约束优化问题，想要求解最优解，可以使用`fmincon`函数。

以下是一个使用`fmincon`函数的示例：```matlabx0=[0,0];%初始点A = []; b = []; Aeq = []; beq = []; % 约束条件lb = [-10, -10]; ub = [10, 10]; % 取值范围options = optimoptions('fmincon'); % 设定优化选项```除了优化选项，Matlab中还有多个参数可供调整，例如算法迭代次数、容差等。

可以根据具体问题的复杂性来调整这些参数。

总而言之，Matlab提供了丰富的最优化工具箱，可以灵活地实现不同类型的最优化方法。

优化方法matlab

优化方法matlab

优化方法matlab对于matlab代码的优化，可以从以下几个方面入手：1. 算法优化：首先，对于算法的优化是最直接有效的方法。

通过优化算法，可以减少代码执行的时间和内存占用。

在编写代码时，可以使用更高效的算法来解决问题。

例如，对于排序问题可以使用快速排序算法代替冒泡排序算法；对于查找问题可以使用二分查找算法代替顺序查找算法。

通过选择合适的算法，可以大大提高程序的效率。

2. 向量化操作：向量化操作是matlab中常用的优化方法之一。

在matlab中，向量和矩阵操作是高效的，而循环操作是低效的。

所以，尽量使用向量和矩阵操作，避免使用循环。

例如，可以使用矩阵乘法代替循环逐个相乘，使用矩阵的元素操作代替循环逐个操作。

3. 减少内存占用：在编写matlab代码时，要注意减少内存的占用，避免不必要的内存拷贝和创建大量的临时变量。

可以使用in-place操作来减少内存使用，尽量避免为临时变量重新分配内存空间。

此外，可以使用matlab内置的函数来高效地处理矩阵和数组，避免不必要的内存开销。

4. 编译优化：matlab提供了mex函数，可以将matlab代码编译成二进制mex 文件，提高代码的执行速度。

通过编译优化，可以将matlab代码转化成C/C++代码，并拥有与C/C++相当的执行效率。

可以将matlab中的瓶颈函数使用mex进行编译优化，提高程序的运行速度。

5. 并行计算：对于一些需要进行大规模计算的问题，可以使用matlab中的并行计算工具箱来进行并行计算，提高程序的运行效率。

可以使用parfor循环来代替普通的for循环，让代码并行执行。

同时，可以使用matlab的并行计算工具箱提供的函数来进行并行计算，如parallel.pool.Constant类来创建共享的常量，parallel.pool.DataQueue类来进行数据通信等。

除了以上几个方面，还可以通过以下方式进行matlab代码的优化：6. 预分配矩阵空间：在编写matlab代码时，可以提前预分配矩阵的空间，避免动态扩展矩阵的大小。

大连理工大学MATLAB概述

大连理工大学MATLAB概述
工作 空间
命令窗口
历史 命令
命令提示符
工作路径
MATLAB变量命名规则
变量名、函数名对字母的大小写是敏感的。 变量名第一个字母必须是英文字母。 变量名可以包含英文字母、下划线和数字。 变量名不能包含空格、标点。
MATLAB预定义变量
预定义变量名 ans eps pi
Inf或inf i或j
数据类型 int8, int16, int32, int64 uint8, uint16, uint32, uint64 single double logical char cell struct funcion_handle
说明 有符号整数 无符号整数 单精度浮点型 双精度浮点型 逻辑型 字符型 单元数组型 结构体型 函数句柄型
MATLAB数值表示
缺省的数据类型为双精度浮点型
例如：3 -10 0.001 1.3e10 1.256e-6 基本操作
ceil( ), floor(), round() %取整
single( )
%单精度浮点型
double( )
%双精度浮点型
MATLAB四则运算符
运算 加 减 乘 除 幂
例：计算sin(45ْ ) >>sin(45*pi/180)
Matalb中正弦函数sin就是常见的正弦函数。 它的参数值是以“弧度”为单位的。 Matlab对字母大小写是敏感的。
例：计算 2ex0.5 1 的值，其中x=4.92。
>>sqrt(2*exp(4.92+0.5)+1)
Matalb中开平方—sqrt(x)，是英文square root的缩写 。 Matalb中指数函数exp(x)，常见的表达方式。

大连理工大学优化方法上机作业

大连理工大学优化方法上机作业

大连理工大学优化方法上机作业本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.March优化方法上机大作业学院：电子信息与电气工程学部姓名：学号：指导老师：上机大作业（一）%目标函数function f=fun(x)f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;end%目标函数梯度function gf=gfun(x)gf=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; End%目标函数Hess矩阵function He=Hess(x)He=[1200*x(1)^2-400*x(2)+2,-400*x(1);-400*x(1), 200;];end%线搜索步长function mk=armijo(xk,dk)beta=0.5; sigma=0.2;m=0; maxm=20;while (m<=maxm)if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk) mk=m; break;endm=m+1;endalpha=beta^mknewxk=xk+alpha*dkfk=fun(xk)newfk=fun(newxk)%最速下降法function [k,x,val]=grad(fun,gfun,x0,epsilon)%功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值maxk=5000; %最大迭代次数beta=0.5; sigma=0.4;k=0;while(k<maxk)gk=feval(gfun,x0); %计算梯度dk=-gk; %计算搜索方向if(norm(gk)<epsilon), break;end%检验终止准则m=0;mk=0;while(m<20) %用Armijo搜索步长if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx0=x0+beta^mk*dk;k=k+1;endx=x0;val=feval(fun,x0);>> x0=[0;0];>> [k,x,val]=grad('fun','gfun',x0,1e-4)迭代次数：k =1033x =0.99990.9998val =1.2390e-008%牛顿法x0=[0;0];ep=1e-4;maxk=10;k=0;while(k<maxk)gk=gfun(x0);if(norm(gk)<ep)x=x0miny=fun(x)k0=kbreak;elseH=inv(Hess(x0));x0=x0-H*gk;k=k+1;endendx =1.00001.0000miny =4.9304e-030迭代次数k0 =2%BFGS方法function [k,x,val]=bfgs(fun,gfun,x0,varargin) %功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值N=1000;epsilon=1e-4;beta=0.55;sigma=0.4;n=length(x0);Bk=eye(n);k=0;while(k<N)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon), break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+beta^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<=oldf+sigma*beta^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+beta^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});>> x0=[0;0];>> [k,x,val]=bfgs('fun','gfun',x0)k =20x =1.00001.0000val =2.2005e-011%共轭梯度法function [k,x,val]=frcg(fun,gfun,x0,epsilon,N)if nargin<5,N=1000;endif nargin<4, epsilon=1e-4;endbeta=0.6;sigma=0.4;n=length(x0);k=0;while(k<N)gk=feval(gfun,x0);itern=k-(n+1)*floor(k/(n+1));itern=itern+1;if(itern==1)dk=-gk;elsebetak=(gk'*gk)/(g0'*g0);dk=-gk+betak*d0; gd=gk'*dk;if(gd>=0),dk=-gk;endendif(norm(gk)<epsilon),break;endm=0;mk=0;while(m<20)if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx=x0+beta^m*dk;g0=gk; d0=dk;x0=x;k=k+1;endval=feval(fun,x);>> x0=[0;0];[k,x,val]=frcg('fun','gfun',x0,1e-4,1000)k =122x =1.00011.0002val =7.2372e-009上机大作业（二）%目标函数function f_x=fun(x)f_x=4*x(1)-x(2)^2-12;%等式约束条件function he=hf(x)he=25-x(1)^2-x(2)^2;end%不等式约束条件function gi_x=gi(x,i)switch icase 1gi_x=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;case 2gi_x=x(1);case 3gi_x=x(2);otherwiseend%求目标函数的梯度function L_grad=grad(x,lambda,cigma)d_f=[4;2*x(2)];d_g(:,1)=[-2*x(1);-2*x(2)];d_g(:,2)=[10-2*x(1);10-2*x(2)];d_g(:,3)=[1;0];d_g(:,4)=[0;1];L_grad=d_f+(lambda(1)+cigma*hf(x))*d_g(:,1);for i=1:3if lambda(i+1)+cigma*gi(x,i)<0L_grad=L_grad+(lambda(i+1)+cigma*gi(x,i))*d_g(:,i+1);continueendend%增广拉格朗日函数function LA=lag(x,lambda,cee)LA=fun(x)+lambda(1)*hf(x)+0.5*cee*hf(x)^2;for i=1:3LA=LA+1/(2*cee)*(min(0,lambda(i+1)+cee*gi(x,i))^2-lambda(i+1)^2); endfunction xk=BFGS(x0,eps,lambda,cigma)gk=grad(x0,lambda,cigma);res_B=norm(gk);k_B=0;a_=1e-4;rho=0.5;c=1e-4;length_x=length(x0);I=eye(length_x);Hk=I;while res_B>eps&&k_B<=10000dk=-Hk*gk;m=0;while m<=5000if lag(x0+a_*rho^m*dk,lambda,cigma)-lag(x0,lambda,cigma)<=c*a_*rho^m*gk'*dkmk=m;break;endm=m+1;endak=a_*rho^mk;xk=x0+ak*dk;delta=xk-x0;y=grad(xk,lambda,cigma)-gk;Hk=(I-(delta*y')/(delta'*y))*Hk*(I-(y*delta')/(delta'*y))+(delta*delta')/(delta'*y);k_B=k_B+1;x0=xk;gk=y+gk;res_B=norm(gk);end%增广拉格朗日法function val_min=ALM(x0,eps)lambda=zeros(4,1);cigma=5;alpha=10;k=1;res=[abs(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);while res>eps&&k<1000xk=BFGS(x0,eps,lambda,cigma);lambda(1)=lambda(1)+cigma*hf(xk);for i=1:3lambda(i+1)=lambda(i+1)+min(0,lambda(i+1)+gi(x0,1)); endk=k+1;cigma=alpha*cigma;x0=xk;res=[norm(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);endval_min=fun(xk);fprintf('k=%d\n',k);fprintf('fmin=%.4f\n',val_min);fprintf('x=[%.4f;%.4f]\n',xk(1),xk(2));>> x0=[0;0];>> val_min=ALM(x0,1e-4)k=10fmin=-31.4003x=[1.0984;4.8779]val_min =-31.4003上机大作业（三）A=[1 1;-1 0;0 -1];n=2;b=[1;0;0];G=[0.5 0;0 2];c=[2 4];cvx_solver sdpt3cvx_beginvariable x(n)minimize (x'*G*x-c*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -2.40.40000.6000A=[2 1 1;1 2 3;2 2 1;-1 0 0;0 -1 0;0 0 -1]; n=3;b=[2;5;6;0;0;0];C=[-3 -1 -3];cvx_solver sdpt3cvx_beginvariable x(n)minimize (C*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -5.40.20000.00001.600011。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

function [x,dk,k]=fjqx(x,s) flag=0;a=0;b=0;k=0;d=1;while(flag==0)[p,q]=getpq(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;x=x+d*s;flag=1;endk=k+1;if(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend%定义求函数值的函数fun，当输入为x0=（x1，x2）时，输出为f function f=fun(x)f=(x(2)-x(1)^2)^2+(1-x(1))^2;function gf=gfun(x)gf=[-4*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),2*(x(2)-x(1)^2)]; function [p,q]=getpq(x,d,s)p=fun(x)-fun(x+d*s)+0.20*d*gfun(x)*s';q=gfun(x+d*s)*s'-0.60*gfun(x)*s';结果：x=[0,1];s=[-1,1];[x,dk,k]=fjqx(x,s)x =-0.0000 1.0000dk =1.1102e-016k =54function f= fun( X )%所求问题目标函数f=X(1)^2-2*X(1)*X(2)+2*X(2)^2+X(3)^2+ X(4)^2-X(2)*X(3)+2*X(1)+3*X(2)-X(3);endfunction g= gfun( X )%所求问题目标函数梯度g=[2*X(1)-2*X(2)+2,-2*X(1)+4*X(2)-X(3)+3,2*X(3)-X(2)-1,2*X(4)];endfunction [ x,val,k ] = frcg( fun,gfun,x0 )%功能：用FR共轭梯度法求无约束问题最小值%输入：x0是初始点，fun和gfun分别是目标函数和梯度%输出：x、val分别是最优点和最优值，k是迭代次数maxk=5000;%最大迭代次数rho=0.5;sigma=0.4;k=0;eps=10e-6;n=length(x0);while(k<maxk)g=feval(gfun,x0);%计算梯度itern=k-(n+1)*floor(k/(n+1));itern=itern+1;%计算搜索方向if(itern==1)d=-g;elsebeta=(g*g')/(g0*g0');d=-g+beta*d0;gd=g'*d;if(gd>=0.0)d=-g;endendif(norm(g)<eps)break;endm=0;mk=0;while(m<20)if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m;break;endm=m+1;endx0=x0+rho^mk*d;val=feval(fun,x0);g0=g;d0=d;k=k+1;endx=x0;val=feval(fun,x0);end结果：>> x0=[0,0,0,0];>> [ x,val,k ] = frcg( 'fun','gfun',x0 )x =-4.0000 -3.0000 -1.0000 0val =-8.0000k =21或者function [x,f,k]=second(x)k=0;dk=dfun(x);g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);while(norm(g1)>=0.02)if(k==3)k=0;g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);else if(k<3)u=((norm(g1))^2)/(norm(g0)^2); s=-g1+u*s;k=k+1;g0=g1;dk=dfun(x);x=x+dk*s;g1=gfun(x);endendf=fun(x);endfunction f=fun(x)f=x(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2+x(4)^2-x(2)*x(3)+2*x(1)+3*x(2)-x(3); function gf=gfun(x)gf=[2*x(1)-2*x(2)+2,-2*x(1)+4*x(2)-x(3)+3,2*x(3)-x(2)-1,2*x(4)];function [p,q]=con(x,d)ss=-gfun(x);p=fun(x)-fun(x+d*ss)+0.2*d*gfun(x)*(ss)';q=gfun(x+d*ss)*(ss)'-0.6*gfun(x)*(ss)';function dk=dfun(x)flag=0;a=0;d=1;while(flag==0)[p,q]=con(x,d);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endEnd结果：x=[0,0,0,0];>> [x,f,k]=second(x)x =-4.0147 -3.0132 -1.0090 0 f = -7.9999k = 1function [f,x,k]=third_1(x)k=0;g=gfun(x);while(norm(g)>=0.001)s=-g;dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)'; j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s)%获取步长flag=0;a=0;d=1;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend结果：x=[0,1];[f,x,k]=third_1(x)f =0.7729x = -0.4196 0.0001k =8（1）程序：function [f,x,k]=third_2(x)k=0;H=inv(ggfun(x));g=gfun(x);while(norm(g)>=0.001)s=(-H*g')';dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2); function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^2);4*x(1)*x(2)*exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2)]; function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)% 步长获取flag=0;a=0;d=1;b=10000;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendEnd结果：x=[0,1];[f,x,k]=third_2(x)f =0.7729x = -0.4193 0.0001k =8（2）程序：function [f,x,k]=third_3(x) k=0;X=cell(2);g=cell(2);X{1}=x;H=eye(2);g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});while(norm(g{2})>=0.001)dx=X{2}-X{1};dg=g{2}-g{1};v=dx/(dx*dg')-(H*dg')'/(dg*H*dg'); h1=H*dg'*dg*H/(dg*H*dg');h2=dx'*dx/(dx*dx');h3=dg*H*dg'*v'*v;H=H-h1+h2+h3;k=k+1;X{1}=X{2};g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});norm(g{2});f=fun(x);x=X{2};endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];function ggf=ggfun(x)ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^2);4*x(1)*x(2)* exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2);function [p,q]=con(x,d,s)p=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';q=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)flag=0;a=0;d=1;b=10000;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0) dk=d;flag=1;endif(p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendend结果：x=[0,1];[f,x,k]=third_3(x)f =0.7729x = -0.4195 0.0000 k=6function callqpactH=[2 0; 0 2];c=[-2 -5]';Ae=[ ]; be=[ ];Ai=[1 -2; -1 -2; -1 2;1 0;0 1];bi=[-2 -6 -2 0 0]';x0=[0 0]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilon=1.0e-9; err=1.0e-6;k=0; x=x0; n=length(x); kmax=1.0e3;ne=length(be); ni=length(bi); lamk=zeros(ne+ni,1); index=ones(ni,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsilon), index(i)=0; endendwhile(k<=kmax)Aee=[];if(ne>0), Aee=Ae; endfor(j=1:ni)if(index(j)>0), Aee=[Aee; Ai(j,:)]; end endgk=H*x+c;[m1,n1] = size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1)); if(norm(dk)<=err)y=0.0;if(length(lamk)>ne)[y,jk]=min(lamk(ne+1:length(lamk))); endif(y>=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk) index(i)=0; break;endendendk=k+1;elseexitflag=1;alpha=1.0; tm=1.0;for(i=1:ni)if((index(i)==0)&(Ai(i,:)*dk<0)) tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk); if(tm1<tm)tm=tm1; ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1), index(ti)=1; end endif(exitflag==0), break; endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x; output.iter=k;function [x,lambda]=qsubp(H,c,Ae,be) ginvH=pinv(H);[m,n]=size(Ae);if(m>0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;end结果>>callqpactx =1.40001.7000lambda =0.8000exitflag =output =fval: -6.4500iter: 7function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)%功能: 用乘子法解一般约束问题: min f(x), s.t. h(x)=0, g(x).=0%输入: x0是初始点, fun, dfun分别是目标函数及其梯度；% hf, dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；% gf, dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；%输出: x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量; % output是结构变量, 输出近似极小值f, 迭代次数, 内迭代次数等maxk=500;c=2.0;eta=2.0;theta=0.8;k=0;ink=0;epsilon=0.00001;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while(btak>epsilon&&k<maxk)%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c);ink=ink+ik;he=feval(hf,x);gi=feval(gf,x);btak=0;for i=1:lbtak=btak+he(i)^2;end%更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp^2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endfor i=1:mlambda(i)=max(0,lambda(i)-c*gi(i));endk=k+1;btaold=btak;x0=x;endendf=feval(fun,x);output.fval=f;output.iter=k;%增广拉格朗日函数function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c) f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for i=1:lpsi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2*c);%不等式约束函数文件g1.mfunction gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;%目标函数的梯度文件df1.mfunction g=df1(x)g=[4, -2*x(2)]';%等式约束（向量）函数的Jacobi矩阵（转置）文件dh1.m function dhe=dh1(x)dhe=[-2*x(1), -2*x(2)]'%不等式约束（向量）函数的Jacobi矩阵（转置）文件dg1.m function dgi=dg1(x)dgi=[10-2*x(1), 10-2*x(2)]';function [x,val,k]=bfgs(fun,gfun,x0,varargin)maxk=500;rho=0.55;sigma=0.4;epsilon=0.00001;k=0;n=length(x0);Bk=eye(n);while(k<maxk)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon)break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});结果x=[2 2]';[x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0) x =1.00134.8987mu =0.7701lambda =0.9434output =fval: -31.9923iter: 4f=[3,1,1];A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0];x=linprog(f,A,b,zeros(3),[0,0,0]',lb)结果：Optimization terminated.x =0.00000.50000.5000。

相关文档
最新文档