大连理工大学优化方法上机大作业程序
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。
教学大纲-大连理工大学教务处
目录《机械设计基础A》 (1)《机械设计基础B》 (8)《**模型设计概论》 (15)阅后删除:请以学部下设学院为单位将全部课程编辑在同一个文档内《机械设计基础A》教学大纲(学分4 学时64)一、课程说明(200字以内,简单说明本课程的地位及教学内容等,阅后删除红色字体)本课程是工科近机械类(包括机械类某些专业)和非机械类专业大类课程之一,是工科学生学习和掌握各种类型的机械中常用机构和通用机械零件的基本知识和基本设计方法的技术基础课。
该课程也是工科学生将来学习专业机械设备课程的理论基础。
本课程在教学内容方面着重基本知识、基本理论和基本设计方法的讲解;在培养实践能力方面着重设计构思和基本设计技能的基本训练。
二、课程目标(对应毕业要求:1-○1、1-○2、1-○3)1. 学习机械工程基础知识和基本理论知识,掌握常用机构的结构、特性等基本知识,了解各种机械的传动原理,具有分析、选用和设计机械设备中基本机构的能力(对应毕业要求:1-○1);2. 通用机械零件的设计原理、方法和机械设计等的一般规律,具有设计机械传动装置和简单机械的能力(对应毕业要求:1-○1);3. 掌握基本的机械设计创新方法,培养学生追求创新的态度和意识(对应毕业要求:1-○1);4. 培养学生树立正确的设计思想,了解机械设计过程中国家有关的经济、环境、法律、安全、健康、伦理等政策和制约因素(对应毕业要求:1-○1);5. 培养学生的工程实践学习能力,使学生掌握典型零件的实验方法,获得实验技能的基本训练,具有运用标准、规范、手册、图册和查阅有关技术资料的能力(对应毕业要求:1-○1);6. 了解机械设计的前沿和新发展动向(对应毕业要求:1-○1)。
三、教学内容、基本要求与学时分配序号教学内容教学要求学时教学方式对应课程目标1 一、基本概念1. 研究的对象、内容;2. 机械设计的基本要求和一般设计过程。
1. 了解本课程研究的对象、内容2. 了解机械设计的基本要求、一般设计过程。
优化方法大作业1
优化方法大作业一、Wolfe-Powell法利用MATLAB软件编写,其中初始值x0=-5;其他参数按照已知条件来取。
当b分别等于8、9、10时,均得到如下结果:而当初值x0变化时,则结果变化比较大,如将x0取-6,则计算结果如下:通过比较可以看出,b的取值对计算结果的影响较小,而初始值x0的取值则对结果影响很大。
从中也表明Wolfe-Powell准则的收敛条件比较弱,容易出现当函数还没取极小值而迭代循环已结束的情况。
具体代码见附录1.二、无约束优化1.DFP法精度ε = 0.02,确定λ使用的是一维非精确搜索算法(直接法,Wolfe-Powell准则)。
结果如下图所示:2. BFGS方法同上一种方法,精度ε = 0. 02,确定λ使用的是一维非精确搜索算法(直接法,Wolfe-Powell准则)。
结果如下图所示:比较两种方法,查看每一步的函数值,得出如下结果。
图1:DFP与BFGS法结果对比图通过图1与表1的比较,BFGS比DFP法最初收敛得更快,但是DFP法比BFGS法的最终结果更好。
DFP与BFGS法的代码分别见附录2、3。
三、Rockafellar乘子法文件myfun.m:function y=myfun(x)y=1000-x(1)^2-2*x(2)^2-x(3)^2-x(1)*x(2)-x(1)*x(3);文件mycon.m:function [c,ceq]=mycon(x)c(1)=(-x(1));c(2)=(-x(2));c(3)=(-x(3));ceq(1)=x(1)^2+x(2)^2+x(3)^2-25;ceq(2)=8*x(1)+14*x(2)+7*x(3)-56文件Q3.m:clearclcx0=[2 2 2]';[x,fval,exitflag,output]=fmincon(@myfun,x0,[],[],[],[],[],[],@mycon);附录1:%% Wolfe-Powell 准则方法clear;tic;c1=0.1;c2=0.65;a=0;b=8;syms x;f=x*x-2*x+7; %函数G0=jacobian(f,x);%求梯度x0=-6; %初始取x=-5% 由于只有一个未知数,默认S=1lambda = 1;x=x0;for k=1:1:1000f0 = eval(f); %求出函数值grad1= eval(G0); %求出梯度值figure1(k) = f0; % 用于绘出迭代过程中函数值变化x= x0+lambda;f1 = eval(f); %求出函数值grad2 = eval(G0); %求出梯度值value1 = f0 - f1 + c1* lambda * grad1 ;value2 = grad2 -c2*grad1 ;if value1 >=-1e-6 && value2 >=-1e-6disp('The variable matrix is : ');disp(x);fun=eval(f);fprintf('The minimum value of function is %f \n',fun);break; %如果满足两个条件,则退出循环。
大连理工大学概率上机作业
大连理工大学概率上机作业第一次上机作业1.利用Matlab自带命令产生1000个均匀随机变量服从U(0,1)。
>> unifrnd(0,1,20,50)ans =Columns 1 through 100.8147 0.6557 0.4387 0.7513 0.3517 0.1622 0.1067 0.8530 0.7803 0.54700.9058 0.0357 0.3816 0.2551 0.8308 0.7943 0.9619 0.6221 0.3897 0.29630.1270 0.8491 0.7655 0.5060 0.5853 0.3112 0.0046 0.3510 0.2417 0.74470.9134 0.9340 0.7952 0.6991 0.5497 0.5285 0.7749 0.5132 0.4039 0.18900.6324 0.6787 0.1869 0.8909 0.9172 0.1656 0.8173 0.4018 0.0965 0.68680.0975 0.7577 0.4898 0.9593 0.2858 0.6020 0.8687 0.0760 0.1320 0.18350.2785 0.7431 0.4456 0.5472 0.7572 0.2630 0.0844 0.2399 0.9421 0.36850.5469 0.3922 0.6463 0.1386 0.7537 0.6541 0.3998 0.1233 0.9561 0.62560.9575 0.6555 0.7094 0.1493 0.3804 0.6892 0.2599 0.1839 0.5752 0.78020.9649 0.1712 0.7547 0.2575 0.5678 0.7482 0.8001 0.2400 0.0598 0.08110.1576 0.7060 0.2760 0.8407 0.0759 0.4505 0.4314 0.4173 0.2348 0.92940.9706 0.0318 0.6797 0.2543 0.0540 0.0838 0.9106 0.0497 0.3532 0.77570.9572 0.2769 0.6551 0.8143 0.5308 0.2290 0.1818 0.9027 0.8212 0.48680.4854 0.0462 0.1626 0.2435 0.7792 0.9133 0.2638 0.9448 0.0154 0.43590.8003 0.0971 0.1190 0.9293 0.9340 0.1524 0.1455 0.4909 0.0430 0.44680.1419 0.8235 0.4984 0.3500 0.1299 0.8258 0.1361 0.4893 0.1690 0.30630.4218 0.6948 0.9597 0.1966 0.5688 0.5383 0.8693 0.3377 0.6491 0.50850.9157 0.3171 0.3404 0.2511 0.4694 0.9961 0.5797 0.9001 0.7317 0.51080.7922 0.9502 0.5853 0.6160 0.0119 0.0782 0.5499 0.3692 0.6477 0.81760.9595 0.0344 0.2238 0.4733 0.3371 0.4427 0.1450 0.1112 0.4509 0.7948Columns 11 through 200.6443 0.3111 0.0855 0.0377 0.0305 0.0596 0.1734 0.9516 0.0326 0.25180.3786 0.9234 0.2625 0.8852 0.7441 0.6820 0.3909 0.9203 0.5612 0.29040.8116 0.4302 0.8010 0.9133 0.5000 0.0424 0.8314 0.0527 0.8819 0.61710.5328 0.1848 0.0292 0.7962 0.4799 0.0714 0.8034 0.7379 0.6692 0.26530.3507 0.9049 0.9289 0.0987 0.9047 0.5216 0.0605 0.2691 0.1904 0.82440.9390 0.9797 0.7303 0.2619 0.6099 0.0967 0.3993 0.42280.3689 0.98270.8759 0.4389 0.4886 0.3354 0.6177 0.8181 0.5269 0.5479 0.4607 0.73020.5502 0.1111 0.5785 0.6797 0.8594 0.8175 0.4168 0.9427 0.9816 0.3439 0.6225 0.2581 0.2373 0.1366 0.8055 0.7224 0.6569 0.4177 0.1564 0.5841 0.5870 0.4087 0.4588 0.7212 0.5767 0.1499 0.6280 0.9831 0.8555 0.1078 0.2077 0.5949 0.9631 0.1068 0.1829 0.6596 0.2920 0.3015 0.6448 0.9063 0.3012 0.2622 0.5468 0.6538 0.2399 0.5186 0.4317 0.7011 0.3763 0.8797 0.4709 0.6028 0.5211 0.4942 0.8865 0.9730 0.0155 0.6663 0.1909 0.8178 0.2305 0.7112 0.2316 0.7791 0.0287 0.6490 0.9841 0.5391 0.4283 0.2607 0.8443 0.2217 0.4889 0.7150 0.4899 0.8003 0.1672 0.6981 0.4820 0.5944 0.1948 0.1174 0.6241 0.9037 0.1679 0.4538 0.1062 0.6665 0.1206 0.0225 0.2259 0.2967 0.6791 0.8909 0.9787 0.4324 0.3724 0.1781 0.5895 0.4253 0.1707 0.3188 0.3955 0.3342 0.7127 0.8253 0.1981 0.1280 0.2262 0.3127 0.2277 0.4242 0.3674 0.6987 0.5005 0.0835 0.4897 0.9991 0.3846 0.1615 0.4357 0.5079 0.9880 0.1978 0.4711 0.1332 0.3395 0.1711 0.5830 0.1788Columns 21 through 300.4229 0.7788 0.2548 0.1759 0.6476 0.5822 0.4046 0.3477 0.8217 0.5144 0.0942 0.4235 0.2240 0.7218 0.6790 0.5407 0.4484 0.1500 0.4299 0.8843 0.5985 0.0908 0.6678 0.4735 0.6358 0.8699 0.3658 0.5861 0.8878 0.5880 0.4709 0.2665 0.8444 0.1527 0.9452 0.2648 0.7635 0.2621 0.3912 0.1548 0.6959 0.1537 0.3445 0.3411 0.2089 0.3181 0.6279 0.0445 0.7691 0.1999 0.6999 0.2810 0.7805 0.6074 0.7093 0.1192 0.7720 0.7549 0.3968 0.4070 0.6385 0.4401 0.6753 0.1917 0.2362 0.9398 0.9329 0.2428 0.8085 0.7487 0.0336 0.5271 0.0067 0.7384 0.1194 0.6456 0.9727 0.4424 0.7551 0.8256 0.0688 0.4574 0.6022 0.2428 0.6073 0.4795 0.1920 0.6878 0.3774 0.7900 0.3196 0.8754 0.3868 0.9174 0.4501 0.6393 0.1389 0.35920.2160 0.3185 0.5309 0.5181 0.9160 0.2691 0.4587 0.5447 0.6963 0.7363 0.7904 0.5341 0.6544 0.9436 0.0012 0.7655 0.6619 0.6473 0.0938 0.3947 0.9493 0.0900 0.4076 0.6377 0.4624 0.1887 0.7703 0.5439 0.5254 0.6834 0.3276 0.1117 0.8200 0.9577 0.4243 0.2875 0.3502 0.7210 0.5303 0.7040 0.6713 0.1363 0.7184 0.2407 0.4609 0.0911 0.6620 0.5225 0.8611 0.4423 0.4386 0.6787 0.9686 0.6761 0.7702 0.5762 0.4162 0.9937 0.4849 0.0196 0.8335 0.4952 0.5313 0.2891 0.3225 0.6834 0.8419 0.2187 0.3935 0.3309 0.7689 0.1897 0.3251 0.6718 0.7847 0.5466 0.8329 0.1058 0.6714 0.4243 0.1673 0.4950 0.1056 0.6951 0.4714 0.4257 0.2564 0.1097 0.7413 0.2703 0.8620 0.1476 0.6110 0.0680 0.0358 0.6444 0.6135 0.0636 0.5201 0.1971 0.9899 0.0550Columns 31 through 400.8507 0.7386 0.5523 0.1239 0.7378 0.5590 0.1781 0.8949 0.6311 0.6925 0.5606 0.5860 0.6299 0.4904 0.0634 0.8541 0.3596 0.0715 0.0899 0.5567 0.9296 0.2467 0.0320 0.8530 0.8604 0.3479 0.0567 0.2425 0.0809 0.3965 0.6967 0.6664 0.6147 0.8739 0.9344 0.4460 0.5219 0.0538 0.7772 0.0616 0.5828 0.0835 0.3624 0.2703 0.9844 0.0542 0.3358 0.4417 0.9051 0.78020.8154 0.6260 0.0495 0.2085 0.8589 0.1771 0.1757 0.0133 0.5338 0.33760.8790 0.6609 0.4896 0.5650 0.7856 0.6628 0.2089 0.8972 0.1092 0.60790.9889 0.7298 0.1925 0.6403 0.5134 0.3308 0.9052 0.1967 0.8258 0.74130.0005 0.8908 0.1231 0.4170 0.1776 0.8985 0.6754 0.0934 0.3381 0.10480.8654 0.9823 0.2055 0.2060 0.3986 0.1182 0.4685 0.3074 0.2940 0.12790.6126 0.7690 0.1465 0.9479 0.1339 0.9884 0.9121 0.45610.7463 0.54950.9900 0.5814 0.1891 0.0821 0.0309 0.5400 0.1040 0.1017 0.0103 0.48520.5277 0.9283 0.0427 0.1057 0.9391 0.7069 0.7455 0.9954 0.0484 0.89050.4795 0.5801 0.6352 0.1420 0.3013 0.9995 0.7363 0.3321 0.6679 0.79900.8013 0.0170 0.2819 0.1665 0.2955 0.2878 0.5619 0.2973 0.6035 0.73430.2278 0.1209 0.5386 0.6210 0.3329 0.4145 0.1842 0.0620 0.5261 0.05130.4981 0.8627 0.6952 0.5737 0.4671 0.4648 0.5972 0.2982 0.7297 0.07290.9009 0.4843 0.4991 0.0521 0.6482 0.7640 0.2999 0.0464 0.7073 0.08850.5747 0.8449 0.5358 0.9312 0.0252 0.8182 0.1341 0.5054 0.7814 0.79840.8452 0.2094 0.4452 0.7287 0.8422 0.1002 0.2126 0.7614 0.2880 0.9430Columns 41 through 500.6837 0.7894 0.1123 0.6733 0.0986 0.9879 0.5975 0.7593 0.8092 0.75190.1321 0.3677 0.7844 0.4296 0.1420 0.1704 0.3353 0.7406 0.7486 0.22870.7227 0.2060 0.2916 0.4517 0.1683 0.2578 0.2992 0.7437 0.1202 0.06420.1104 0.0867 0.6035 0.6099 0.1962 0.3968 0.4526 0.1059 0.5250 0.76730.1175 0.7719 0.9644 0.0594 0.3175 0.0740 0.4226 0.6816 0.3258 0.67120.6407 0.2057 0.4325 0.3158 0.3164 0.6841 0.3596 0.4633 0.5464 0.71520.3288 0.3883 0.6948 0.7727 0.2176 0.4024 0.5583 0.2122 0.3989 0.64210.6538 0.5518 0.7581 0.6964 0.2510 0.9828 0.7425 0.0985 0.4151 0.41900.7491 0.2290 0.4326 0.1253 0.8929 0.4022 0.4243 0.8236 0.1807 0.39080.5832 0.6419 0.6555 0.1302 0.7032 0.6207 0.4294 0.1750 0.2554 0.81610.7400 0.4845 0.1098 0.0924 0.5557 0.1544 0.1249 0.1636 0.0205 0.31740.2348 0.1518 0.9338 0.0078 0.1844 0.3813 0.0244 0.6660 0.9237 0.81450.7350 0.7819 0.1875 0.4231 0.2120 0.1611 0.2902 0.8944 0.6537 0.78910.9706 0.1006 0.2662 0.6556 0.0773 0.7581 0.3175 0.5166 0.9326 0.85230.8669 0.2941 0.7978 0.7229 0.9138 0.8711 0.6537 0.7027 0.1635 0.50560.0862 0.2374 0.4876 0.5312 0.7067 0.3508 0.9569 0.1536 0.9211 0.63570.3664 0.5309 0.7690 0.1088 0.5578 0.6855 0.9357 0.9535 0.7947 0.95090.3692 0.0915 0.3960 0.6318 0.3134 0.2941 0.4579 0.5409 0.5774 0.44400.6850 0.4053 0.2729 0.1265 0.1662 0.5306 0.2405 0.6797 0.4400 0.06000.5979 0.1048 0.0372 0.1343 0.6225 0.8324 0.7639 0.0366 0.2576 0.8667 2.参考课本综合例题2.5.4和2.5.5中的方法,模拟产生1000个随机变量,使其服从参数为2的指数分布,进而计算这1000个随机数的均值和方差。
大连理工优化方法大作业MATLAB编程
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。
大连理工大学c++大作业
大连理工大学C++程序设计大作业班级:111111111姓名:暗暗暗暗暗暗暗暗学号:1111111111邮箱:1111111111111111任课教师:赵国辉上交时间:2013.7.22目录1.第一次上机作业 (3)2.第二次上机作业 (6)3.第三次上机作业 (10)4.第四次上机作业 (15)5.结课大作业 (18)6.课堂总结 (48)1. 第一次上机作业1.1作业要求1.整型和浮点型的二进制表示2.一个整数的二进制(64位)转换成十进制表示。
1.2 核心代码说明整型和浮点型的二进制表示#include <iostream>#include <cstdio>#include <cstring>#include <algorithm>#include <queue>#include <stack>using namespace std;void getint(int x){stack<int>s;while(x){s.push(x&1);x>>=1;}while(!s.empty()){cout<<s.top();s.pop();}} //将整型转换成二进制函数void getfloat(float y){queue<int>q;int x=(int)y;getint(x);y-=x; //将浮点型转换成二进制函数if(!y) return ;putchar('.');while(y){y*=2.0;if(y>=1.0){q.push(1);y-=1.0;}else q.push(0);}while(!q.empty()){cout<<q.front();q.pop();}}int main(){int x;float y;cin>>x;//输入一个整数getint(x);cout<<endl;cin>>y;//输入一个浮点数getfloat(y);cout<<endl;return 0;}一个整数的二进制(64位)转换成十进制表示。
大连理工大学优化方法上机作业
大连理工大学优化方法上机作业本页仅作为文档页封面,使用时可以删除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。
大连理工优化方法大作业MATLAB编程
fun ctio n [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)A2)A2+(1-x(1)F2;function gf=gfun(x)gf=[-4*x(1)*(x (2) -x(1)A2)+2*(x(1)-1),2*(x(2)-x(1)A2)];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 =54取初始= (0.0. 0,0)r^'l用兵柜梯皮法求解下面无约東优化问题:min f (x) = x孑—2x^X2 十2x孑 + x孑H-爲—X2天3 十 2xj + 3|X2 —*3,其中步长g的选取可利用习題1戎精确一维披索.注:通过比习题验证共範梯度法求辉门无二次西数极小点至多需要“次迭代.fun ctio n f= fun( X )%所求问题目标函数f=X(1)A2-2*X(1)*X (2)+2*X(2)A2+X(3)A2+ X(4) A2-X( 2)*X(3)+2*X(1)+3*X(2)-X(3);end function 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)];end function [ x,val,k ] = frcg( fun,gfun,xO )%功能:用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(fu n,xO+rhoAm*d)<feval(fu n,xO)+sigma*rhoAm*g'*d) mk=m; break ;endm=m+1;endx0=x0+rho A 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 =或者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))A2)/( norm(gO)A2); 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(1F2-2*x(1)*x (2)+2*x (2)A2+x(3)A2+x(4)A2-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 = 1取初始点3 = (0」)二考虑下面无约東优化问题:min f(x)二冷 + 2x2 + exp(xf + 天孑),其中歩长Qk的选取可別用习题1或精确一维搜索•搜索方向为一HNW ♦取垃=b•取皿=R2f防)]"9耳丈啟为BFG5公式亠通过此习题体会上述三种算法的收敛速度.fun ctio n [f,x,k]=third_1(x) k=0;g=gfu n(x);while (norm(g)>=0.001) s=-g;dk=dfu n( x,s);x=x+dk*s;k=k+1;g=gfu n(x);f=fun( x);endfun ctio n f=fun(x)f=x(1)+2*x(2)A2+exp(x(1)A2+x(2)A2);fun ctio n gf=gfu n(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)]; 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}; end结果: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)A2+exp(x(1F2+x(2)A2);function gf=gfun(x) gf=[1+2*x(1)*exp(x(1F2+x(2)A2),4*x(2)+2*x(2)*(x(1F2+x(2)A2)]; function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x (2) A2),4*x(1)*x (2) *exp(x(1)A2+x(2)A2);4*x(1)*x(2)*exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2)];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;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)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)A2+exp(x(1F2+x(2)A2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)];function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x(2)A2),4*x(1)*x(2)*exp(x(1)A2+x(2)A2);4*x(1)*x(2)* exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2);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;if (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.41950.0000 k=6*U 用有效集法求解下面勺勺二次规划问题:(XI 一 I)2 + (x 2 一 2.5)2 X1 - 2X2 + 2 > 0-Xi — 2>(2 + 6 > 0-Xi + 2X2 + 2 > 0xi,x 2 > 0function 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,xO)fun ctio n [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilo n=1.0e-9; err=1.0e-6;k=0; x=x0; n=len gth(x); kmax=1.0e3;n e=le ngth(be); ni=le ngth(bi); lamk=zeros( ne+n i,1); in dex=ones(n i,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsil on), i ndex(i)=0; end while (k<=kmax)mmSi.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 ;endendk=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;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(fu n, hf,gf,dfu n, dhf,dgf,xO)%功能:用乘子法解一般约束问题: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;i nk=0;epsilo n=0.00001;x=xO;he=feval(hf,x);gi=feval(gf,x);n=len gth(x);l=le ngth(he);m=le ngth(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(y2;end% 更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp A2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endlambda(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);s仁 s1+he(y2;psi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3A2-lambda(i)A2;endpsi=psi+s2/(2*c);% 不等式约束函数文件 g1.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-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)n ewf=feval(fu n, x0+rho A m*dk,vararg in {:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rhoAm*gk'*dk) mk=m;break ;endm=m+1;endx=x0+rhoAmk*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: 4利用序列二次规划方法求解习题5中的约束优化问题:min 4xi 一好一 12s.t. 25 - x? —x孑=Q10x一召 + 10旳-xj - 34 > 0 X1,X2 > 0tf=[3,1,1];A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0]; x=li nprog(f,A,b,zeros(3),[0,0,0]',lb)结果:Optimization terminated.0.00000.50000.5000。
大连理工大学机械设计大作业
目录一、设计任务书及原始数据 (2)二、根据已知条件计算传动件的作用力 (3)2.1计算齿轮处转矩T、圆周力F t 、径向力F r及轴向力F a .. 3 2.2计算链轮作用在轴上的压力 (3)2.3计算支座反力 (4)三、初选轴的材料,确定材料的机械性能 (4)四、进行轴的结构设计 (5)4.1确定最小直径 (5)4.2设计其余各轴段的直径和长度,且初选轴承型号 (5)4.3选择连接形式与设计细部结构 (6)五.轴的疲劳强度校核 (6)5.1轴的受力图 (6)5.2绘制弯矩图 (7)5.3绘制转矩图 (8)5.4确定危险截面 (9)5.5计算当量应力,校核轴的疲劳强度 (9)六、选择轴承型号,计算轴承寿命 (10)6.1计算轴承所受支反力 (10)6.2计算轴承寿命 (11)七、键连接的计算 (11)八、轴系部件的结构装配图 (12)一、设计任务书及原始数据题目二:二级展开式斜齿圆柱齿轮减速器输出轴组合结构设计表1 设计方案及原始数据二、根据已知条件计算传动件的作用力2.1计算齿轮处转矩T、圆周力F t、径向力F r及轴向力F a已知:轴输入功率P=4.3kW,转速n=130r/(min)。
(1)齿轮上的力转矩计算公式:T=9.550×106P/n将数据代入公式中,得:T=315885(N·mm)圆周力计算公式:Ft=2T/,==416(mm) (认为是法面模数)将转矩T带入其中,得:Ft=1519(N)径向力计算公式:Fr =Ft×tanα/cos,=将圆周力Ft 带入其中,得:Fr=558(N)轴向力计算公式:Fa = Ft×tan将圆周力Ft 带入其中,得:Fa=216(N)2.2计算链轮作用在轴上的压力链轮的分度园直径链速v=链的圆周力F=链轮作用在轴上的压力2.3计算支座反力1、计算垂直面(XOZ)支反力根据受力分析图,我们可以利用垂直面力矩平衡原理(ΣMY=0)得出求解A点垂直面支反力Rz1:R z1= Ft1- Rz2ΣM a= R Z2× l AC- F t1× l AB=0即ΣMa =1519 ×135-RZ2× 215=0RZ2=954 NRz1=565 N2、计算水平面(XOY)支反力根据受力分析图,我们可以利用水平面力矩平衡原理(ΣMZ=0)得出求解A点水平面支反力Ry1的计算公式:R y1= FQ- Ry2–FrΣM Z= R y2× l AC+ F r× l AB +F a×r- F Q× l ADΣM Z=R y2× 215+558×135+216×135-×315=0Ry2=4437NRy1=- 4437 –558=-1635 N三、初选轴的材料,确定材料的机械性能初选材料及机械性能见表四、进行轴的结构设计4.1确定最小直径按照扭转强度条件计算轴的最小值d min。
大连理工大学 秋季优化方法大作业
m=m+1; end x0=x0+rho^mk*d; val=feval(fun,x0); g0=g; d0=d; k=k+1; end x=x0; val=feval(fun,x);
//f(x)//
function f=fun(x) x1=[1 0]*x; x2=[0 1]*x; f=(1-x1)^2+100*(x2-x1^2)^2; //梯度函数// function g=gfun(x) x1=[1 0]*x; x2=[0 1]*x; g=[-2*(1-x1)-400*x1*(x2-x1^2); 200*(x2-x1^2)]; //运行过程// >> x0=[0 0]'
k=k+1; btaold=btak; x0=x; end f=feval(fun,x); output.fval=f; output.iter=k; output.inner_iter=ink; output.bta=btak;
//增广拉格朗日函数//
function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x); l=length(he); m=length(gi); psi=f; s1=0.0; for(i=1:l) psi=psi-he(i)*mu(i); s1=s1+he(i)^2; end psi=psi+0.5*sigma*s1; s2=0.0; for(i=1:m) s3=max(0.0, lambda(i) - sigma*gi(i)); s2=s2+s3^2-lambda(i)^2; end psi=psi+s2/(2.0*sigma); //增广拉格朗日函数// function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) dpsi=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) dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i); end for(i=1:m) dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i); end //f(x)// function f=f1(x) f=4*x(1)-x(2)^2-12; //等式约束// function he=h1(x) he=25-x(1)^2-x(2)^2;
大连理工优化方法大作业MATLAB编程
fun ctio n [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)A2)A2+(1-x(1)F2;function gf=gfun(x)gf=[-4*x(1)*(x (2) -x(1)A2)+2*(x(1)-1),2*(x(2)-x(1)A2)];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 =54取初始= (0.0. 0,0)r^'l用兵柜梯皮法求解下面无约東优化问题:min f (x) = x孑—2x^X2 十2x孑 + x孑H-爲—X2天3 十 2xj + 3|X2 —*3,其中步长g的选取可利用习題1戎精确一维披索.注:通过比习题验证共範梯度法求辉门无二次西数极小点至多需要“次迭代.fun ctio n f= fun( X )%所求问题目标函数f=X(1)A2-2*X(1)*X (2)+2*X(2)A2+X(3)A2+ X(4) A2-X( 2)*X(3)+2*X(1)+3*X(2)-X(3);end function 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)];end function [ x,val,k ] = frcg( fun,gfun,xO )%功能:用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(fu n,xO+rhoAm*d)<feval(fu n,xO)+sigma*rhoAm*g'*d) mk=m; break ;endm=m+1;endx0=x0+rho A 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 =或者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))A2)/( norm(gO)A2); 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(1F2-2*x(1)*x (2)+2*x (2)A2+x(3)A2+x(4)A2-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 = 1取初始点3 = (0」)二考虑下面无约東优化问题:min f(x)二冷 + 2x2 + exp(xf + 天孑),其中歩长Qk的选取可別用习题1或精确一维搜索•搜索方向为一HNW ♦取垃=b•取皿=R2f防)]"9耳丈啟为BFG5公式亠通过此习题体会上述三种算法的收敛速度.fun ctio n [f,x,k]=third_1(x) k=0;g=gfu n(x);while (norm(g)>=0.001) s=-g;dk=dfu n( x,s);x=x+dk*s;k=k+1;g=gfu n(x);f=fun( x);endfun ctio n f=fun(x)f=x(1)+2*x(2)A2+exp(x(1)A2+x(2)A2);fun ctio n gf=gfu n(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)]; 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}; end结果: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)A2+exp(x(1F2+x(2)A2);function gf=gfun(x) gf=[1+2*x(1)*exp(x(1F2+x(2)A2),4*x(2)+2*x(2)*(x(1F2+x(2)A2)]; function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x (2) A2),4*x(1)*x (2) *exp(x(1)A2+x(2)A2);4*x(1)*x(2)*exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2)];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;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)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)A2+exp(x(1F2+x(2)A2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)];function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x(2)A2),4*x(1)*x(2)*exp(x(1)A2+x(2)A2);4*x(1)*x(2)* exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2);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;if (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.41950.0000 k=6*U 用有效集法求解下面勺勺二次规划问题:(XI 一 I)2 + (x 2 一 2.5)2 X1 - 2X2 + 2 > 0-Xi — 2>(2 + 6 > 0-Xi + 2X2 + 2 > 0xi,x 2 > 0function 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,xO)fun ctio n [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilo n=1.0e-9; err=1.0e-6;k=0; x=x0; n=len gth(x); kmax=1.0e3;n e=le ngth(be); ni=le ngth(bi); lamk=zeros( ne+n i,1); in dex=ones(n i,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsil on), i ndex(i)=0; end while (k<=kmax)mmSi.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 ;endendk=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;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(fu n, hf,gf,dfu n, dhf,dgf,xO)%功能:用乘子法解一般约束问题: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;i nk=0;epsilo n=0.00001;x=xO;he=feval(hf,x);gi=feval(gf,x);n=len gth(x);l=le ngth(he);m=le ngth(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(y2;end% 更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp A2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endlambda(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);s仁 s1+he(y2;psi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3A2-lambda(i)A2;endpsi=psi+s2/(2*c);% 不等式约束函数文件 g1.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-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)n ewf=feval(fu n, x0+rho A m*dk,vararg in {:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rhoAm*gk'*dk) mk=m;break ;endm=m+1;endx=x0+rhoAmk*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: 4利用序列二次规划方法求解习题5中的约束优化问题:min 4xi 一好一 12s.t. 25 - x? —x孑=Q10x一召 + 10旳-xj - 34 > 0 X1,X2 > 0tf=[3,1,1];A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0]; x=li nprog(f,A,b,zeros(3),[0,0,0]',lb)结果:Optimization terminated.0.00000.50000.5000。
优化作业流程的三种方法
优化作业流程的三种方法一、简化步骤。
1.1 去除冗余操作。
在作业流程里啊,常常有些步骤就像“画蛇添足”一样,没什么实际用处,还拖慢速度。
比如说,有些文件填写,要求填一大堆重复信息,这不是浪费时间嘛。
咱就得把这些多余的部分找出来,大刀阔斧地砍掉。
就像整理房间,把那些没用的杂物都扔掉,房间立马就清爽了。
1.2 合并相似任务。
有些任务啊,就像双胞胎似的,很相似。
那咱们就别分开做了,把它们合并起来。
例如,数据收集和初步整理,这俩事儿关联性很强,完全可以一起做。
这样就像把散落在各处的小珠子串成一条项链,一下子就整齐有序了,效率也能大大提高。
二、引入新技术。
2.1 自动化工具的使用。
现在科技这么发达,有好多自动化的工具就像“得力助手”一样。
像一些办公软件,能自动处理文档格式、进行数据计算。
以前人工做这些得花老长时间了,现在点几下鼠标就搞定。
比如说财务报表制作,有了专门的软件,数据输入进去,各种图表、分析立马就出来了,又快又准,简直是“神来之笔”。
2.2 信息化管理系统。
建立信息化管理系统那可是个好办法。
这就好比给作业流程装上了一个智能大脑。
所有的任务分配、进度跟踪、资源调配都能在这个系统里一目了然。
员工们也能清楚地知道自己该干什么,什么时候干。
就像大家都在一条清晰的轨道上行驶的列车,不会乱套。
三、员工培训与激励。
3.1 针对性培训。
员工要是对作业流程不熟悉,就像盲人摸象一样,只能瞎干。
所以要进行针对性的培训。
根据不同岗位的需求,把流程详细地教给员工。
这就像给战士们配上合适的武器并且教会他们怎么用。
比如生产线上的工人,要让他们清楚每个环节的操作规范,这样才能保证产品质量,提高生产效率。
3.2 激励机制。
人都是需要激励的,有了激励就像给汽车加足了油。
设立一些奖励制度,对于那些在作业流程优化方面有好点子或者执行得特别好的员工,给予奖励。
可以是奖金、荣誉称号或者晋升机会。
这样员工们就会像打了鸡血一样,积极主动地去寻找优化作业流程的方法,整个团队也就充满了活力。
优化方法课程大作业
优化方法课程上机大作业学部:电子信息与电气工程学部专业:生物医学工程班级:电信硕1303学号:21309210姓名:史益新大连理工大学Dalian University of Technology解:(1)MATLAB代码如下:clc;clear all;close all;[x,y]=meshgrid(-2:0.1:2,-1:0.1:3);z=(y-x.^2).^2+(1-x).^2;mesh(x,y,z)hold on;xk=[0 1]';epsilon=1e-5;plot(xk(1),xk(2),'ro');text(xk(1),xk(2),'start point');hold on;[ x,val,k ]=Newton('fun','gfun','Hess',xk,epsilon)plot(x(1),x(2),'ro');text(x(1),x(2),'end point');function [ x,val,k ] = Newton( fun,gfun,Hess,xk,epsilon ) k=0;while(1)gk=feval(gfun,xk);hk=feval(Hess,xk);sk=-inv(hk)*gk;if(norm(gk)<epsilon)break;endxk=xk+sk;k=k+1;endx=xk;val=feval(fun,x);endfunction f = fun(x)f=(x(2)-x(1)^2)^2+(1-x(1))^2;endfunction g=gfun(x)g=[-4*x(1)*(x(2)-x(1)^2)+2*(1-x(1)),2*(x(2)-x(1)^2)]'; endfunction He = Hess( x )n=length(x);He=zeros(n,n);He=[12*x(1)^2-4*x(2)+2,-4*x(1);-4*x(1), 2 ];end(2)代码运行结果如下:解:(1)MATLAB代码如下:clc;clear all;close all;x0=[0 0 0 0]';epsilon=1e-5;[ x,val,k ]=Frcg('fun','gfun',x0,epsilon)function [ x,val,k ] = Frcg( fun,gfun,x0,epsilon )rho=0.6;sigma=0.5;k=0;n=length(x0);while(1)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)<epsilon)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,x);function 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)]';end(2)代码运行结果如下:解:(1)MATLAB代码如下:clc;clear all;close all;[x,y]=meshgrid(-2:0.1:2,-1:0.1:3);z=5*(y-x.^2).^2+(x-1).^2;mesh(x,y,z)hold on;x0=[2 0]';epsilon=1e-5;plot(x0(1),x0(2),'ro');text(x0(1),x0(2),'start point');hold on;[ x1,val1,k1 ]=grad('fun','gfun',x0,epsilon)plot(x1(1),x1(2),'ro');text(x1(1),x1(2),'end point1');hold on;[ x2,val2,k2 ]=znNewton('fun','gfun','Hess',x0,epsilon) plot(x2(1),x2(2),'bo');text(x2(1),x2(2),'end point2');hold on;[ x3,val3,k3 ]=BFGS('fun','gfun',x0,epsilon)plot(x3(1),x3(2),'go');text(x3(1),x3(2),'end point3');hold on;function [ x,val,k ] = grad( fun,gfun,x0,epsilon )rho=0.5;sigma=0.4;k=0;while(1)g=feval(gfun,x0);d=-g;if(norm(d)<epsilon)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;k=k+1;endx=x0;val=feval(fun,x0);endfunction [ x,val,k ] = znNewton( fun,gfun,Hess,x0,epsilon )rho=0.5;sigma=0.4;k=0;while(1)gk=feval(gfun,x0);hk=feval(Hess,x0);dk=-inv(hk)*gk;if(norm(gk)<epsilon)break;endm=0;mk=0;while(m<20)if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk) mk=m;break;endm=m+1;endx0=x0+rho^mk*dk;k=k+1;endx=x0;val=feval(fun,x0);endfunction [ x,val,k ] = BFGS( fun,gfun,x0,epsilon )rho=0.5;sigma=0.4;k=0;n=length(x0);bk=eye(n);while(1)gk=feval(gfun,x0);if(norm(gk)<epsilon)break;enddk=-inv(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);endk=k+1;x0=x;endval=feval(fun,x0);endfunction f = fun(x)f=5*(x(2)-x(1)^2)^2+(x(1)-1)^2;endfunction g=gfun(x)g=[-20*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),10*(x(2)-x(1)^2)]'; endfunction He = Hess( x )n=length(x);He=zeros(n,n);He=[60*x(1)^2-20*x(2)+2,-20*x(1);-20*x(1), 10 ];end(2)代码运行结果如下:解:(1)MATLAB程序如下:clc;close all;clear all;x0=[1,0]';epsilon=1e-5;[ x,mu,lambda,output ] = multphr( 'f1','h1','g1','df1','dh1','dg1',x0,epsilon )function [ x,mu,lambda,output ] = multphr( fun,hf,gf,dfun,dhf,dgf,x0,epsilon )%MULTPHR Summary of this function goes here% Detailed explanation goes heresigma=2.0;eta=2.0;theta=0.8;k=0;ink=0;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(gi);%initial of multi-vectormu=0.1*ones(l,1);lambda=0.1*ones(m,1);btak=10;btaold=10;while(btak>epsilon)[x,ival,ik]=BFGS('mpsi','dmpsi',x0,epsilon,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma);ink=ink+ik;he=feval(hf,x);gi=feval(gf,x);btak=0;for(i=1:l)btak=btak+he(i)^2;endfor(i=1:m)temp=min(gi(i),lambda(i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);if(btak>epsilon)if(k>=2 & btak>theta*btaold)sigma=eta*sigma;endfor(i=1:l)mu(i)=mu(i)-sigma*he(i);endfor(i=1:m)lambda(i)=max(0,lambda(i)-sigma*gi(i));endendk=k+1;btaold=btak;x0=x;endf=feval(fun,x);output.fval=f;output.iter=k;output.inner_iter=ink;output.bta=btak;endfunction [ x,val,k ] = BFGS( fun,gfun,x0,varargin )rho=0.5;epsilon=1e-5;sigma=0.4;k=0;n=length(x0);bk=eye(n);while(1)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon)break;enddk=-inv(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{:});endfunction psi = mpsi( x,epsilon,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma ) %MPSI Summary of this function goes here% Detailed explanation goes heref=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for(i=1:l)psi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*sigma*s1;s2=0;for(i=1:m)s3=max(0,lambda(i)-sigma*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2*sigma);endfunction dpsi = dmpsi( x,epsilon,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma ) %DMPSI Summary of this function goes here% Detailed explanation goes heredpsi=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)dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);endfor(i=1:m)dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);endendfunction f = f1( x )f=4*x(1)-x(2)^2-12;endfunction gi = g1( x )gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;%unequation constrainendfunction he = h1( x )he=25-x(1)^2-x(2)^2;%equation constrainendfunction g = df1( x )g=[4,-2*x(2)]';endfunction dgi = dg1( x )dgi=[-2*x(1)+10,-2*x(2)+10]';endfunction dhe = dh1( x )dhe=[-2*x(1),-2*x(2)]';end(2)代码运行结果如下:解:(1)MATLAB代码如下:clc;close all;clear all;H=[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]';epsilon=1e-9;[ x,lamk,exitflag,output ] = qpact( H,c,Ae,be,Ai,bi,x0,epsilon )function [ x,lamk,exitflag,output ] = qpact( H,c,Ae,be,Ai,bi,x0,epsilon ) %QPACT Summary of this function goes here% Detailed explanation goes hereerr=1e-6;k=0;x=x0;n=length(x);kmax=1e3;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,:)];endendgk=H*x+c;[m1,n1]=size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1));if(norm(dk)<=err)y=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;tm=1;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;end %updata the setk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;endfunction [ x,lambda ] = qsubp( H,c,Ae,be )%QSUBP Summary of this function goes here % Detailed explanation goes hereginvH=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;endend(2)代码运行结果如下:解:(1)MATLAB代码如下:clc;close all;clear all;x0=[0.5 0.2]';mu0=[ ]';lam0=[0 0 0 0]';epsilon=1e-6;[ x,mu,lam,val,k ] = sqpm( x0,mu0,lam0,epsilon)function [ x,mu,lam,val,k ] = sqpm( x0,mu0,lam0,epsilon ) %SQPM Summary of this function goes here% Detailed explanation goes heren=length(x0);l=length(mu0);m=length(lam0);rho=0.5;eta=0.1;B0=eye(n);x=x0;mu=mu0;lam=lam0;Bk=B0;sigma=0.8;[hk,gk]=cons(x);dfk=df1(x);[Ae,Ai]=dcons(x);Ak=[Ae;Ai];k=0;while(1)[dk,mu,lam]=qpsubp(dfk,Bk,Ae,hk,Ai,gk,epsilon);mp1=norm(hk,1)+norm(max(-gk,0),1);if (norm(dk,1)<epsilon) & (mp1<1e-5)break;enddeta=0.05;tau=max(norm(mu,inf),norm(lam,inf));if(sigma*(tau+deta)<1)sigma=sigma;elsesigma=1.0/(tau+2*deta);endim=0;while(im<=20)if(phi1(x+rho^im*dk,sigma)-phi1(x,sigma)<eta*rho^im*dphi1(x,sigma,dk)) mk=im;break;endim=im+1;if(im==20)mk=10;endendalpha=rho^mk;x1=x+alpha*dk;[hk,gk]=cons(x1);dfk=df1(x1);[Ae,Ai]=dcons(x1);Ak=[Ae;Ai];lamu=pinv(Ak)'*dfk;if(l>0 & m>0)mu=lamu(1:l);lam=lamu(l+1:l+m);endif(l==0)mu=[];lam=lamu;endif(m==0)mu=lamu;lam=[];endsk=alpha*dk;yk=dlax(x1,mu,lam)-dlax(x,mu,lam);if(sk'*yk>0.2*sk'*Bk*sk)theta=1;elsetheta=0.8*sk'*Bk*sk/(sk'*Bk*sk-sk'*yk);endzk=theta*yk+(1-theta)*Bk*sk;Bk=Bk+zk*zk'/(sk'*zk)-(Bk*sk)*(Bk*sk)'/(sk'*Bk*sk);x=x1;k=k+1;endval=f1(x);endfunction [ d,mu,lam,val,k ] = qpsubp( dfk,Bk,Ae,hk,Ai,gk,epsilon ) %QPSUBP Summary of this function goes here% Detailed explanation goes heren=length(dfk);l=length(hk);m=length(gk);gamma=0.05;rho=0.5;sigma=0.2;ep0=0.05;mu0=0.05*zeros(l,1);lam0=0.05*zeros(m,1);d0=ones(n,1);u0=[ep0;zeros(n+l+m,1)];z0=[ep0;d0;mu0;lam0,];k=0;z=z0;ep=ep0;d=d0;mu=mu0;lam=lam0;while(1)dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);if(norm(dh)<epsilon)break;endA=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);b=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)*u0-dh;dz=pinv(A)*b;if(l>0 & m>0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);dl=dz(n+l+2:n+l+m+1);endif(l==0)de=dz(1);dd=dz(2:n+1);dl=dz(n+2:n+m+1);endif(m==0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);endi=0;while(i<=20)if(l>0 & m>0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(l==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(m==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam,dfk,Bk,Ae,hk,Ai,gk);endif(norm(dh1)<=(1-sigma*(1-gamma*ep0)*rho^i)*norm(dh))mk=i;break;endi=i+1;if(i==20)mk=10;endendalpha=rho^mk;if(l>0 & m>0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;lam=lam+alpha*dl;endif(l==0)ep=ep+alpha*de;d=d+alpha*dd;lam=lam+alpha*dl;endif(m==0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;endk=k+1;endendfunction dh = dah( ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk )%DAH Summary of this function goes here% Detailed explanation goes heren=length(dfk);l=length(hk);m=length(gk);dh=zeros(n+l+m+1,1);dh(1)=ep;if(l>0 & m>0)dh(2:n+1)=Bk*d-Ae'*mu-Ai'*lam+dfk;dh(n+2:n+l+1)=hk+Ae*d;for(i=1:m)dh(n+l+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);endendif(l==0)dh(2:n+1)=Bk*d-Ai'*lam+dfk;for(i=1:m)dh(n+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);endendif(m==0)dh(2:n+1)=Bk*d-Ae'*mu+dfk;dh(n+2:n+l+1)=hk+Ae*d;enddh=dh(:);endfunction bet = beta( ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma )%BETA Summary of this function goes here% Detailed explanation goes heredh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);bet=gamma*norm(dh)*min(1,norm(dh));endfunction [ h,g ] = cons( x )%CONS Summary of this function goes here% Detailed explanation goes hereh=[ ];g=[-x(1)^2+6*x(1)-4*x(2)+11,x(1)*x(2)-3*x(2)-exp(x(1)-1)+1,x(1),x(2)]'; endfunction [ dh,dg ] = dcons( x )%DCONS Summary of this function goes here% Detailed explanation goes heredh=[ ];dg=[-2*x(1)+6,-4;x(2)-exp(x(1)-1),x(1)-1;1,0;0,1];endfunction [ dd1,dd2,v1 ] = ddv( ep,d,lam,Ai,gk)%DDV Summary of this function goes here% Detailed explanation goes herem=length(gk);dd1=zeros(m,m);dd2=zeros(m,m);v1=zeros(m,1);for(i=1:m)fm=sqrt(lam(i)^2+(gk(i)+Ai(i,:)*d)^2+2*ep^2);dd1(i,i)=1-lam(i)/fm;dd2(i,i)=1-(gk(i)+Ai(i,:)*d)/fm;v1(i)=-2*ep/fm;endendfunction df = df1( x )%DF1 Summary of this function goes here% Detailed explanation goes heredf=[2*x(1)-16,2*x(2)-10]';endfunction dl = dlax( x,mu,lam )%DLAX Summary of this function goes here% Detailed explanation goes heredf=df1(x);[Ae,Ai]=dcons(x);[m1,m2]=size(Ai);[l1,l2]=size(Ae);if(l1==0)dl=df-Ai'*lam;endif(m1==0)dl=df-Ae'*mu;endif(l1>0 & m1>0)dl=df-Ae'*mu-Ai'*lam;endendfunction f = f1( x )%F1 Summary of this function goes here% Detailed explanation goes heref=x(1)^2+x(2)^2-16*x(1)-10*x(2);endfunction p = phi( ep,a,b )%PHI Summary of this function goes here% Detailed explanation goes herep=a+b-sqrt(a^2+b^2+2*ep^2);endfunction p = phi1( x,sigma )%PHI1 Summary of this function goes here% Detailed explanation goes heref=f1(x);[h,g]=cons(x);gn=max(-g,0);l0=length(h);m0=length(g);if(l0==0)p=f+1.0/sigma*norm(gn,1);endif(m0==0)p=f+1.0/sigma*norm(h,1);endif(l0>0 & m0>0)p=f+1.0/sigma*(norm(h,1)+norm(gn,1)); endendfunction dp = dphi1( x,sigma,d )%DPHI1 Summary of this function goes here% Detailed explanation goes heredf=df1(x);[h,g]=cons(x);gn=max(-g,0);l0=length(h);m0=length(g);if(l0==0)dp=df'*d-1.0/sigma*norm(gn,1);endif(m0==0)dp=df'*d-1.0/sigma*norm(h,1);endif(l0>0 & m0>0)dp=df'*d-1.0/sigma*(norm(h,1)+norm(gn,1)); endendfunction A = JacobiH( ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk )%JACOBIH Summary of this function goes here% Detailed explanation goes heren=length(dfk);l=length(hk);m=length(gk);A=zeros(n+l+m+1,n+l+m+1);[dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk);if(l>0 & m>0)A=[1, zeros(1,n), zeros(1,l), zeros(1,m);zeros(n,1), Bk, -Ae', -Ai';zeros(l,1), Ae, zeros(l,l), zeros(l,m);v1, dd2*Ai, zeros(m,l), dd1]; endif(l==0)A=[1, zeros(1,n), zeros(1,m);zeros(n,1), Bk, -Ai';v1, dd2*Ai, dd1];endif(m==0)A=[1, zeros(1,n), zeros(1,l);zeros(n,1), Bk, -Ae';zeros(l,1), Ae, zeros(l,l)];endend(2)代码运行结果如下:解:(1)MATLAB代码如下:clc;close all;clear all;f=[1 1 1 1 1 1 1 1];A=[1 0 0 0.5;0 1 0.2 0.3;0 0.1 1 0.2];Aeq=[A,-A];beq=[-1 0.2 1]';lb=zeros(8,1);x=linprog(f,[],[],Aeq,beq,lb);x=[x(1)-x(5) x(2)-x(6) x(3)-x(7) x(4)-x(8)]' (2)代码运行结果如下:。
大连理工大学庞丽萍最优化方法MATLAB程序
班级:优化1班授课老师:庞丽萍姓名:学号:第二章12.(1)用修正单纯形法求解下列LP问题:>>clear>>A=[121100;123010;215001];[m,n]=size(A);b=[10;15;20];r=[-1-2-31];c=[-1-2-31];bs=[3:3];nbs=[1:4];a1=A(:,3);T=A(:,bs);a2=inv(T)*a1;b=inv(T)*b;A=[eye(m),a2];B=eye(m);xb=B\b;cb=c(bs);cn=c(nbs);con=1;M=zeros(1);while conM=M+1;t=cb/B;r=c-t*A;if all(r>=0)x(bs)=xb;x(nbs)=0;fx=cb*xb;disp(['当前解是最优解,minz=',num2str(fx)])disp('对应的最优解为,x=')disp(x)breakendrnbs=r(nbs);kk=find(rnbs==min(rnbs));k=kk(1);Anbs=A(:,nbs);yik=B\Anbs(:,k);xb=B\b;%yi0if all(yik<=0)disp('此LP问题无有限的最优解,计算结束',x)disp(xb)breakelsei=find(yik>0);w=abs(xb(i,1)./yik(i,1));l=find(w==min(w));rr=min(l);yrrk=yik(rr,1);Abs=A(:,bs);D=Anbs(:,k);Anbs(:,k)=Abs(:,rr);Abs(:,rr)=D;F=bs(rr);bs(rr)=nbs(k);nbs(k)=F;AA=[Anbs,Abs];EE=eye(m);EE(:,rr)=-yik./yrrk;Errk=EE;Errk(rr,rr)=1/yrrk;BB=Errk/B;B=inv(BB);cb=c(:,bs);xb=Errk*xb;x(bs)=xb;x(nbs)=0;fx=cb*xb;endif M>=1000disp('此问题无有限最优解')breakendend%结果当前解是最优解,minz=-15对应的最优解为,x=2.5000 2.5000 2.50000第三章30题DFP算法求函数极小点的计算程序function[x,val,k]=dfp(fun,gfun,x0)%功能:用DFP算法求解无约束问题:minf(x)%输入:x0是初始点,fun,gfun分别是目标函数及其梯度%输出:x,val分别是近似最优点和最优值,k是迭代次数.maxk=1e5;%给出最大迭代次数rho=0.55;sigma=0.4;epsilon=1e-5;k=0;n=length(x0);Hk=inv(feval('Hess',x0));%Hk=eye(n);while(k<maxk)gk=feval(gfun,x0);%计算梯度if(norm(gk)<epsilon),break;end%检验终止准则dk=-Hk*gk;%解方程组,计算搜索方向m=0;mk=0;while(m<20)%用Armijo搜索求步长if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk’*dk)mk=m;break;endm=m+1;end%DFP校正x=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x)-gk;if(sk'*yk>0)Hk=Hk-(Hk*yk*yk'*Hk)/(yk'*Hk*yk)+(sk*sk')/(sk'*yk);endk=k+1;x0=x;endval=feval(fun,x0);%习题26的程序调用方式及结果:function y=fun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=(x(1)-1)^2+5*(x2-x(1)^2)^2endfunction y=gfun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[20]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=1.000001.00000val=k=6%习题27的程序调用方式及结果:function y=fun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=x1+2*x(2)^2+exp(x(1)^2+x(2)^2)endfunction y=gfun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[10]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=-0.419360val=0.77291k=536题编写Hooke-Jeeves方法求函数极小点的计算程序。
优化方法MATLAB编程——大连理工大学
The optimal solution is 0.000000.
The optimal "x" is "ans".
ans =
1.0000 1.0000 1.0000 1.0000 可以看出,用Newton法经过6次迭代就能求出最优解。
3. BFGS法 源程序如下: function zy_x=di2tiBFGS(x) %该函数用来解大作业第二题。 x0=x; yimuxulong=1e-5; k=0; g0=g(x0); H0=eye(4);s0=-H0*g0;
3
s0=-g0;
大连理工大学优化方法上机大作业
end end x=x0+lanmed*s0; x0=x; g0=g(x); s0=-g0; k=k+1; end end zy_x=x; zyj=f(x); fprintf('after %d iterations,obtain the optimal solution.\n\nThe optimal solution is %f.\n\n The optimal "x" is "ans".\n',k,zyj);
4
大连理工大学优化方法上机大作业
>> x=[-1.2 1 -1.2 1]'; >> di2titidu(x) after 5945 iterations,obtain the optimal solution.
The optimal solution is 0.000000.
The optimal "x" is "ans".
大连理工大学优化方法上机大作业
大连理工大学《人工智能》大作业及要求
学习中间:专业:年级:学号:学生:题目:1.谈谈你对本课程学习过程中的心得当会与主张?经过这门课程的学习,我对人工智能有了一些简略的理性知道,我晓得了人工智能从诞生到开展阅历一个绵长的过程,许多人为此做出了不懈的尽力。
我觉得这门课程是一门赋有应战性的科学,而从事这项工作的人不只要懂得计算机常识,还需求懂得心思学和哲学。
2. 《人工智能》课程设计, 从以下5个题目中任选其一作答。
《人工智能》课程设计留意:从以下5个题目中任选其一作答。
总则:不约束编程语言,提交word文档,不要提交紧缩包作业提交:大作业上交时文件名写法为:[名字奥鹏卡号学习中间](如:戴卫东101410013979浙江台州奥鹏学习中间[1]VIP)以附件word文档方式上交离线作业(附件的巨细约束在10M以内),挑选已完结的作业(留意命名),点提交即可。
如下图所示。
留意事项:独立完结作业,禁绝抄袭其别人或许请人代做,如有相同作业,分数以零分计!题目一:A*算法要求:(1)编撰一份word文档,里边包含(算法思路、算法程序框图、重排九宫疑问)章节。
(2)算法思路:简略介绍该算法的根本思想,100字摆布即可。
(3)算法程序框图:制作流程图或原理图,从算法的开端到完毕的程序框图。
(4)关于重排九宫疑问的启示式函数: f (x)= p(x)+3s(x)p(x)是x结点和方针结点比较每个将牌“离家”的最短间隔之和;s(x)是:每个将牌和方针比较,若该将牌的后继和方针中该将牌的后继不一样,则该将牌得2分,一样则该将牌得0分,中心方位有将牌得1分,没将牌得0分。
关于给定的初始格式和方针状况请按此启示式函数给出查找的状况空间图。
初始格式方针状况题目二:回归算法要求:(1)编撰一份word文档,里边包含(常见的回归算法、根据实例的算法详细细节)章节。
(2)常见的回归算法包含:最小二乘法(Ordinary Least Square),逻辑回归(Logistic Regression),逐渐式回归(Stepwise Regression),多元自习惯回归样条(Multivariate Adaptive Regression Splines)以及本地散点滑润估量(Locally Estimated Scatterplot Smoothing),请挑选一个算法描绘下算法中心思想(3)随意选用一个实例完成你所挑选的回归算法。
大连理工程序设计 第4次上机作业
1.(1)题目描述:排序:使用起泡法对10个整数实现递减排序。
题目分析:将两个相邻数比较,大的调到前头。
采用两重循环。
有10个数,则要进行9趟比较,第i趟比较中要进行10-i次两两比较。
流程图:程序代码:#include<stdio.h>main(){int a[10],i,j,k,t;printf("please input ten numbers:\n");for(i=0;i<10;i++)scanf("%d",&a[i]);for(k=1;k<10;k++)for(j=0;j<10-k;j++)if(a[j]<a[j+1]){t=a[j];a[j]=a[j+1];a[j+1]=t;}printf("The sorted numbers:\n");for(i=0;i<10;i++)printf("%d ",a[i]);}程序截屏:1.(2)题目描述:排序:使用选择法对10个整数实现递减排序。
题目分析:以冒泡排序法为基础,在两两比较后并不马上进行交换,而是在找到最大的数之后,记住最大数的位置(数组中的下标),待一轮比较完毕后,再将最大的数一次交换到位。
编程时使用二重循环。
外循环控制进行N-1轮排序,内循环找出第i 轮的最大值。
流程图:程序代码:#include<stdio.h>main(){int a[10],i,j,k,m,t;printf("please input ten numbers:\n");for(i=0;i<10;i++)scanf("%d",&a[i]);for(k=0;k<9;k++){t=k;for(j=k+1;j<10;j++)if(a[t]<a[j])t=j;{m=a[k];a[k]=a[t];a[t]=m;}}printf("The sorted numbers:\n");for(i=0;i<10;i++)printf("%d ",a[i]);}程序截屏:2.题目描述:统计:从键盘输入10个学生的数学(MT)、英语(EN)和物理(PH)成绩,并按照如下统计形式输出,包括学生学号(NO)、各科成绩、总成绩(SUM)、平均分(AVE)及是否每科都超过90分(‘Y’or ‘N’)题目分析:利用二维数组,输入十个学生的成绩。
大连理工大学优化方法上机大作业
学院:专业:班级:学号:姓名:上机大作业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 + *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--The 2-th iter, the residual isx1 =法: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 + *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*(fa 0)')/((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<=1000 ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + *ak*slope ak = 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+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=; %下降标准系数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 : 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 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime -------------------------------------------------------------------0|||+01|+00|+02|+01 +00| 0:0:00| chol 1 11|||||+01|+00 +01| 0:0:01| chol 1 12|||||+00|+00 +01| 0:0:01| chol 1 13|||||+00|+00 +00| 0:0:01| chol 1 14||||||+00 +00| 0:0:01| chol 1 15||||||+00 +00| 0:0:01| chol 1 16||||||+00 +00| 0:0:01| chol 1 17||||||+00 +00| 0:0:01| chol 1 18||||||+00 +00| 0:0:01|stop: max(relative gap, infeasibilities) <------------------------------------------------------------------- number of iterations = 8primal objective value = +00dual objective value = +00gap := trace(XZ) =relative gap =actual relative gap =rel. primal infeas (scaled problem) =rel. dual " " " =rel. primal infeas (unscaled problem) = +00rel. dual " " " = +00norm(X), norm(y), norm(Z) = +00, +00, +00norm(A), norm(b), norm(C) = +00, +00, +00Total CPU time (secs) =CPU time per iteration =termination code = 0DIMACS: +00 +00-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval):2、>> clear alln=2; c=[-2,-4]'; G=[,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 : 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 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime -------------------------------------------------------------------0||||+00|+02| +01 +00| 0:0:00| chol 1 11|||||+01| +00 | 0:0:00| chol 1 12|||||+00| +00 | 0:0:00| chol 1 13|||||| | 0:0:00| chol 1 14|||||| | 0:0:00| chol 1 15|||||| | 0:0:00| chol 1 16|||||| | 0:0:00| chol 1 17|||||| | 0:0:00| chol 1 18|||||| | 0:0:00| chol 1 19|||||| | 0:0:00| chol 1 110|||||| | 0:0:00| chol 2 211|||||| | 0:0:00| chol 2 212|||||| | 0:0:00| chol 2 213|||||| | 0:0:00| chol 2 214|||||| | 0:0:00|stop: max(relative gap, infeasibilities) <------------------------------------------------------------------- number of iterations = 14primal objective value =dual objective value =gap := trace(XZ) =relative gap =actual relative gap =rel. primal infeas (scaled problem) =rel. dual " " " =rel. primal infeas (unscaled problem) = +00rel. dual " " " = +00norm(X), norm(y), norm(Z) = +00, +00, +00norm(A), norm(b), norm(C) = +00, +00, +00Total CPU time (secs) =CPU time per iteration =termination code = 0DIMACS: +00 +00-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval): -3。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
函数定义:% 目标函数function f = fun(x)fm=0;for i=1:499fmi = (1-x(i))^2 + 100*(x(i+1)-x(i)^2)^2;fm=fm+fmi;endf =fm;end% 梯度function g = grad(x)g = zeros(500,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2));for i=2:499g(i)=2*(x(i)-1)+400*x(i)*(x(i)^2-x(i+1))+200*(x(i)-x(i-1)^2); endg(500) = 200*(x(500)-x(499)^2);end% 二阶梯度function g = grad2(x)g = zeros(500,500);g(1,1)=2+400*(3*x(1)^2-x(2));g(1,2)=-400*x(1);for i=3:500g(1,i)=0;endfor i=1:498g(500,i)=0;endg(500,499)=-400*x(499);g(500,500)=200;for i=2:499for j=1:500if j==i-1g(i,j)= -400*x(i-1);elseif j==ig(i,j)= 2+400*(3*x(i)^2-x(i+1))+200;elseif j==i+1g(i,j)= -400*x(i);else g(i,j)=0;endendendend1.最速下降法function x_star = steepest(x0,eps)gk = grad(x0);res = norm(gk);k = 0;while res > eps && k<=50000dk = -gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;gk = grad(xk);res = norm(gk);fprintf('--After the %d-th iter, the residual is %f\n',k,res);endx_star = xk;end运行结果:>> steepest(zeros(500,1),0.0001)--After the 1-th iter, the residual is 447.213595--After the 2-th iter, the residual is 164.162920……………………………--After the 13587-th iter, the residual is 0.000111--After the 13588-th iter, the residual is 0.000098ans =1.0000 (以下略,x各解均为1.0000)2.牛顿法(不使用一维线搜索法,由于初始点离结果较远导致不收敛,故采用初始值x均为0.8)function 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('--After the %d-th iter, the residual is %f\n',k,res); endx_star = xk;end运行结果:>> newton(0.8*ones(500,1),0.0001)--After the 1-th iter, the residual is 95613.530371--After the 2-th iter, the residual is 27958.861084--After the 3-th iter, the residual is 8027.532086--After the 4-th iter, the residual is 2197.358772--After the 5-th iter, the residual is 524.525307--After the 6-th iter, the residual is 82.202566--After the 7-th iter, the residual is 4.168551--After the 8-th iter, the residual is 0.010203--After the 9-th iter, the residual is 0.002320--After the 10-th iter, the residual is 0.000000ans =1.0000(以下略,x各解均为1.0000)3.BFGS法function x_star = bfgs(x0,eps)g0 = grad(x0);gk=g0;res = norm(gk);Hk=eye(500);k = 0;while res > eps && k<=10000dk = -Hk*gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x=xk-x0;x0 = xk;g0=gk;gk = grad(xk);g=gk-g0;H0=Hk;Hk=H0+((1+g'*H0*g/(x'*g))*x*x'-H0*g*x'-x*g'*H0)/(x'*g); res = norm(gk);fprintf('--After the %d-th iter, the residual is %f\n',k,res); endx_star = xk;end运行结果:>> bfgs(zeros(500,1),0.0001)--After the 1-th iter, the residual is 447.213595--After the 2-th iter, the residual is 408.218892 ……………………………--After the 321-th iter, the residual is 0.000105--After the 322-th iter, the residual is 0.000076ans =1.0000(以下略,x各解均为1.0000)4.共轭梯度法function x_star =fr(x0,eps)gk = grad(x0);res = norm(gk);k = 0;dk = -gk;while res > eps && k<=10000ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk); while f1 > f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endd0=dk;g0=gk;k=k+1;x0=xk;gk=grad(xk);f=(norm(gk)/norm(g0))^2;res=norm(gk);dk=-gk+f*d0;fprintf('--After the %d-th iter, the residual is %f\n',k,res); endx_star = xk;end运行结果:>> fr(zeros(500,1),0.0001)--After the 1-th iter, the residual is 447.213595--After the 2-th iter, the residual is 434.581110 ………………………………--After the 4288-th iter, the residual is 0.000100--After the 4289-th iter, the residual is 0.000100ans =1.0000(以下略,x各解均为1.0000)1.linprog求解f = [-3;-1;-3];A = [2 1 1;1 2 3;2 2 1];b = [2;5;6];lb = [0;0;0];[x,fval,exitflag,output,lambda] = linprog(f,A,b,[],[],lb) 运行结果:>> aOptimization terminated.x = 0.2000 0.0000 1.6000fval = -5.4000exitflag = 1output =iterations: 7algorithm: 'interior-point'cgiterations: 0message: 'Optimization terminated.'constrviolation: 0firstorderopt: 1.1410e-12lambda =ineqlin: [3x1 double]eqlin: [0x1 double]upper: [3x1 double]lower: [3x1 double]2.CVX求解n=3;a=[-3 -1 -3];b=[2;5;6];C=[2 1 1;1 2 3;2 2 1];cvx_beginvariable x(n)minimize( a*x)subject toC * x <= bx>=0cvx_endx运算结果:Calling 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 1 1|0.912|1.000|9.4e-01|4.6e-02|6.5e+01|-5.606627e+00 -2.967567e+01| 0:0:00| chol 1 12|1.000|1.000|1.3e-07|4.6e-03|8.5e+00|-2.723981e+00 -1.113509e+01| 0:0:00| chol 1 1 3|1.000|0.961|2.3e-08|6.2e-04|1.8e+00|-4.348354e+00 -6.122853e+00| 0:0:00| chol 1 1 4|0.881|1.000|2.2e-08|4.6e-05|3.7e-01|-5.255152e+00 -5.622375e+00| 0:0:00| chol 1 1 5|0.995|0.962|1.6e-09|6.2e-06|1.5e-02|-5.394782e+00 -5.409213e+00| 0:0:00| chol 1 1 6|0.989|0.989|2.7e-10|5.2e-07|1.7e-04|-5.399940e+00 -5.400100e+00| 0:0:00| chol 1 1 7|0.989|0.989|5.3e-11|5.8e-09|1.8e-06|-5.399999e+00 -5.400001e+00| 0:0:00| chol 1 1 8|1.000|0.994|2.8e-13|4.3e-11|2.7e-08|-5.400000e+00 -5.400000e+00| 0:0:00|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.36CPU time per iteration = 0.03termination 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.4x =0.20000.00001.6000。