优化matlab作业,现代设计优化算法MATLAB实现

开篇语

前阵子做现代设计方法的时候&＃xff0c;发现网上很是缺乏这种作业形式的简易算法实现&＃xff0c;所以特地来简书写一篇。有两份&＃xff0c;一份是我的(说来惭愧&＃xff0c;我的大部分都是在网上找的代码&＃xff0c;然后在自己的电脑上跑一次&＃xff0c;跑出来了就行了的。而且我的电脑跑到2-11就扑街了。暂时还没有拿去修&＃xff0c;所以&＃xff0c;其实我的代码都是在网上整理之后调整了下就上的&＃xff0c;准确性不敢保证)另一份是我的室友的&＃xff0c;据他说是全部经过调试的&＃xff0c;虽然还是有不少的错误&＃xff0c;但是应该比我的要好一点。

我的电脑

另外&＃xff0c;我的电脑因为崩了&＃xff0c;所以我的代码无法验证结果&＃xff0c;为了交作业&＃xff0c;只能把我的室友的那些运行结果直接上一遍了。估计有点出入&＃xff0c;见谅&＃xff0c;重要的是代码

正文

2-10

2-10

我的&＃xff1a;

黄金分割法&＃xff1a;

f&＃61;&＃64;(x) x&＃43;20/x

golden(f,2,10,0.01)

function[xmin]&＃61;golden(f,a,b,e)

k&＃61;0;

a1&＃61;b-0.618*(b-a); %插入点的值

a2&＃61;a&＃43;0.618*(b-a);

while b-a>e %循环条件

y1&＃61;subs(f,a1);

y2&＃61;subs(f,a2);

if y1>y2 %比较插入点的函数值的大小

a&＃61;a1; %进行换名

a1&＃61;a2;

y1&＃61;y2;

a2&＃61;a&＃43;0.618*(b-a);

else

b&＃61;a2;

a2&＃61;a1;

y2&＃61;y1;

a1&＃61;b-0.618*(b-a);

end

k&＃61;k&＃43;1;

end %迭代到满足条件为止就停止迭代

xmin&＃61;(a&＃43;b)/2;

fmin&＃61;subs(f,xmin) %输出函数的最优值

fprintf(&＃39;k&＃61;\n&＃39;); %输出迭代次数

disp(k);

f &＃61; &＃64;(x)x&＃43;20/x

>> [x,y]&＃61;golden(f,2,10,0.01)

x &＃61;

4.4683

y &＃61;

8.9443

二次插值法&＃xff1a;

f&＃61;&＃64;(x) x&＃43;20/x;

a&＃61;2;b&＃61;10;

eps&＃61;1.0e-6; % 计算精度

x1&＃61;a;x3&＃61;b;

x2&＃61;(a&＃43;b)/2;

f1 &＃61; f(x1);f2 &＃61; f(x2);f3 &＃61; f(x3);

while 1

C1&＃61;(x2^2-x3^2)*f1&＃43;(x3^2-x1^2)*f2&＃43;(x1^2-x2^2)*f3;

C2&＃61;(x2-x3)*f1&＃43;(x3-x1)*f2&＃43;(x1-x2)*f3;

xp&＃61;0.5*C1/C2;

fp&＃61; f(xp);

if abs(x2-xp)<&＃61;eps % 区间长度小于eps时

if abs(f2-fp)<&＃61;eps % df小于eps时退出

if fp<&＃61;f2

xmin &＃61; xp

fmin &＃61; f(xp) % 极小值

break;

else

xmin &＃61; x2

fmin &＃61; f(x2) % 极小值

break;

end

end

else

if fp<&＃61;f2

if xp<&＃61;x2

x3&＃61;x2;

x2&＃61;xp;

f3&＃61;f2;

f2&＃61;fp;

else

x1&＃61;x2;

x2&＃61;xp;

f1&＃61;f2;

f2&＃61;fp;

end

else

if xp<&＃61;x2

x1&＃61;xp;

f1&＃61;fp;

else

x3&＃61;xp;

f3&＃61;fp;

end

end

end

end

f&＃61;x&＃43;20/x;[xmin,fmin]&＃61;main(f,2,10,0.01)

xmin &＃61;

4.4869

fmin &＃61;

8.9443

室友的&＃xff1a;

黄金分割法函数&＃xff1a;

f&＃61;&＃64;(x) x&＃43;20/x

yellowking(f,2,10,0.01)

function[xmin]&＃61;yellowking(f,a,b,e)

k&＃61;0;

a1&＃61;b-0.618*(b-a); %插入点的值

a2&＃61;a&＃43;0.618*(b-a);

while b-a>e %循环条件

y1&＃61;subs(f,a1);

y2&＃61;subs(f,a2);

if y1>y2 %比较插入点的函数值的大小

a&＃61;a1; %进行换名

a1&＃61;a2;

y1&＃61;y2;

a2&＃61;a&＃43;0.618*(b-a);

else

b&＃61;a2;

a2&＃61;a1;

y2&＃61;y1;

a1&＃61;b-0.618*(b-a);

end

k&＃61;k&＃43;1;

end %迭代到满足条件为止就停止迭代

xmin&＃61;(a&＃43;b)/2;

fmin&＃61;subs(f,xmin) %输出函数的最优值

fprintf(&＃39;k&＃61;\n&＃39;); %输出迭代次数

disp(k);

结果指令&＃xff1a;

f&＃61;&＃64;(x)x&＃43;20/x

f &＃61;

&＃64;(x)x&＃43;20/x

>> [x,y]&＃61;gold(f,2,10,0.01)

x &＃61;

4.4683

y &＃61;

8.9443

二次插值法函数&＃xff1a;

function [xmin,fmin]&＃61; main(f,a0,b0,epsilon)

a&＃61;a0;

b&＃61;b0;

x1&＃61;a;

f1&＃61;f(x1);

x3&＃61;b;

f3&＃61;f(x3);

x2&＃61;5;

f2&＃61;f(x2);

c1&＃61;(f3-f1)/(x3-x1);

c2&＃61;((f2-f1)/(x2-x1)-c1)/(x2-x3);

xp&＃61;0.4*(x1&＃43;x3-c1/c2);fp&＃61;f(xp);

while (abs(xp-x2)>&＃61;epsilon)

if x2

if f2>fp

f1&＃61;f2;x1&＃61;x2;

x2&＃61;xp;f2&＃61;fp;

else

f3&＃61;fp;x3&＃61;xp;

end

else

if f2>fp

f3&＃61;f2;x3&＃61;x2;

f2&＃61;fp;x2&＃61;xp;

else

f1&＃61;fp;x2&＃61;xp;

end

end

c1&＃61;(f3-f1)/(x3-x1);

c2&＃61;((f2-f1)/(x2-x1)-c1)/(x2-x3);

xp&＃61;0.5*(x1&＃43;x3-c1/c2);

fp&＃61;f(xp);

end

if f2>fp

xmin&＃61;xp;fmin&＃61;f(xp);

else

xmin&＃61;x2;fmin&＃61;f(x2);

end

end

结果&＃xff1a;

clear all;f&＃61;x&＃43;20/x;[xmin,fmin]&＃61;main(f,2,10,0.01)

xmin &＃61;

4.4869

fmin &＃61;

8.9443

2-11

2-11

我的&＃xff1a;

2-11

function [k ender]&＃61;steepest(f,x,e)

%梯度下降法,f为目标函数(两变量x1和x2)&＃xff0c;x为初始点,如[3;4]

syms x1 x2 m; %m为学习率

d&＃61;-[diff(f,x1);diff(f,x2)]; %分别求x1和x2的偏导数&＃xff0c;即下降的方向

flag&＃61;1; %循环标志

k&＃61;0; %迭代次数

while(flag)

d_temp&＃61;subs(d,x1,x(1)); %将起始点代入&＃xff0c;求得当次下降x1梯度值

d_temp&＃61;subs(d_temp,x2,x(2)); %将起始点代入&＃xff0c;求得当次下降x2梯度值

nor&＃61;norm(d_temp); %范数

if(nor>&＃61;e)

x_temp&＃61;x&＃43;m*d_temp; %改变初始点x的值

f_temp&＃61;subs(f,x1,x_temp(1)); %将改变后的x1和x2代入目标函数

f_temp&＃61;subs(f_temp,x2,x_temp(2));

h&＃61;diff(f_temp,m); %对m求导&＃xff0c;找出最佳学习率

m_temp&＃61;solve(h); %求方程&＃xff0c;得到当次m

x&＃61;x&＃43;m_temp*d_temp; %更新起始点x

k&＃61;k&＃43;1;

else

flag&＃61;0;

end

end

ender&＃61;double(x); %终点

end

syms x1 x2;

f&＃61;x1^2&＃43;x2^2-x1*x2-10*x1-4*x2&＃43;60;

x&＃61;[0;0];

e&＃61;0.01;

[k ender]&＃61;steepest(f,x,e)

ender &＃61;

7.9961

5.9971

室友的&＃xff1a;

2-11:

梯度函数&＃xff1a;

function [k,ender]&＃61;tidu(f,x,e)

syms x1 x2 m;

d&＃61;-[diff(f,x1);diff(f,x2)];

flag&＃61;1;

k&＃61;0;

while(flag)

d_temp&＃61;subs(d,x1,x(1));

d_temp&＃61;subs(d_temp,x2,x(2));

nor&＃61;norm(d_temp);

if(nor>&＃61;e)

x_temp&＃61;x&＃43;m*d_temp;

f_temp&＃61;subs(f,x1,x_temp(1));

f_temp&＃61;subs(f_temp,x2,x_temp(2));

h&＃61;diff(f_temp,m);

m_temp&＃61;solve(h);

x&＃61;x&＃43;m_temp*d_temp;

k&＃61;k&＃43;1;

else

flag&＃61;0;

end

end

ender&＃61;double(x);

end

结果指令&＃xff1a;

syms x1 x2;

f&＃61;x1^2&＃43;x2^2-x1*x2-10*x1-4*x2&＃43;60;

x&＃61;[0;0];

e&＃61;0.01;

[k ender]&＃61;tidu(f,x,e)

ender &＃61;

7.9961

5.9971

2-12

我的&＃xff1a;

2-12

展开为二阶泰勒式

syms x1 x2;

taylor(x1^4&＃43;2*x2^3-3*x1^2*x2)

ans &＃61;

3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) &＃43; 3*(x1 - 1)^2 &＃43; 6*(x2 - 1)^2 - 1

牛顿法求解&＃xff1a;

function all&＃61;newton(f,x,e)

syms x1 x2 h;

d&＃61;-[diff(f,x1);diff(f,x2)];

h&＃61;hessian(f);

flag&＃61;1;

h1&＃61;h^-1;

while (flag)

d_temp&＃61;subs(d,x1,x(1));

d_temp&＃61;subs(d_temp,x2,x(2));

nor&＃61;norm(d_temp);

if(nor>&＃61;e)

x&＃61;x&＃43;h1*d_temp;

else

flag&＃61;0;

end

end

all&＃61;double(x);

结果指令&＃xff1a;

clear all

>> syms x1 x2;

f&＃61;3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) &＃43; 3*(x1 - 1)^2 &＃43; 6*(x2 - 1)^2 - 1;

x&＃61;[1;1];

e&＃61;0.01;all&＃61;newton(f,x,e)

all &＃61;

1.1667

0.8333

室友的&＃xff1a;

2-12:

展开为二阶泰勒式

syms x1 x2;

taylor(x1^4&＃43;2*x2^3-3*x1^2*x2)

ans &＃61;

3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) &＃43; 3*(x1 - 1)^2 &＃43; 6*(x2 - 1)^2 - 1

牛顿函数&＃xff1a;

function all&＃61;newton(f,x,e)

syms x1 x2 h;

d&＃61;-[diff(f,x1);diff(f,x2)];

h&＃61;hessian(f);

flag&＃61;1;

h1&＃61;h^-1;

while (flag)

d_temp&＃61;subs(d,x1,x(1));

d_temp&＃61;subs(d_temp,x2,x(2));

nor&＃61;norm(d_temp);

if(nor>&＃61;e)

x&＃61;x&＃43;h1*d_temp;

else

flag&＃61;0;

end

end

all&＃61;double(x);

结果指令&＃xff1a;

clear all

>> syms x1 x2;

f&＃61;3*x2 - 2*x1 - 6*(x1 - 1)*(x2 - 1) &＃43; 3*(x1 - 1)^2 &＃43; 6*(x2 - 1)^2 - 1;

x&＃61;[1;1];

e&＃61;0.01;all&＃61;newton(f,x,e)

all &＃61;

1.1667

0.8333

2-13(1)

我的&＃xff1a;

2-13(1)

外点惩罚函数法&＃xff1a;

function [ x,y ] &＃61; Epfm_min( fx,gx,hx,xx0,s,c,a)

%fx是目标函数

%gx是不等式约束方程组(且g>&＃61;0)

%xx0是初始点

%hx是等式约束方程组(且h&＃61;0)

%s是精确度(s>0)

%c是放大系数(c>1)

%a是罚因子(默认为1)

syms x1 x2

xx1&＃61;xx0;

v&＃61;[x1,x2];

a1&＃61;a;

Pxk&＃61;1;%假设Px等于1&＃xff0c;以免不必要错误

G&＃61;-subs(gx,v,xx1);%用于判别max{0,-g(x)}

while Pxk>s

if(G<0)

Px&＃61;a1*hx*hx;

else

Px&＃61;a1*hx*hx&＃43;a1*gx*gx;

end

Fx&＃61;fx&＃43;a1*Px;%将约束问题化为了一个无约束的问题

% 接下来解min F(x)

dFx1&＃61;diff(Fx,x1);%分别对x1&＃xff0c;x2求偏导数

dFx2&＃61;diff(Fx,x2);

[k,b]&＃61;solve(dFx1,dFx2,&＃39;x1&＃39;,&＃39;x2&＃39;);%求出

xx2&＃61;xx1&＃43;[k,b];

Pxk&＃61;a1*subs(Px,v,xx2);

xx1&＃61;xx2;%相当于置k&＃61;k&＃43;1

a1&＃61;c*a1;%罚因子放大

G&＃61;-subs(gx,v,xx1);%用于判别max{0,-g(x)}

end

x&＃61;xx1;

y&＃61;a1/c;

syms x1 x2;

fx&＃61;x1&＃43;x2;

gx&＃61;-x1;

hx&＃61;x1^2-x2;

s&＃61;10.^-5

c&＃61;10

xx0&＃61;[0,0]

a&＃61;1;

[x,y]&＃61;Epfm_min( fx,gx,hx,xx0,s,c,a)

>> x&＃61;[0.1;0.2];k&＃61;0.1;e&＃61;0.01;r&＃61;1;[x,minf]&＃61; Epfm_min (p,x,k,r,e)

x &＃61;

0.0015

minf &＃61;

0.0030

Ans&＃61;

0.0045

内点惩罚函数法

function [x,minf]&＃61;minNF(f,x0,g,u,v,var,eps)

format long;

if nargin&＃61;&＃61;6

eps&＃61;1.0e-4;

end

k&＃61;0;

FE&＃61;0;

for i&＃61;1:length(g)

FE&＃61;FE&＃43;1/g(i);

end

x1&＃61;transpose(x0);

x2&＃61;inf;

while 1

FF&＃61;u*FE;

SumF&＃61;f&＃43;FF;

[x2,minf]&＃61;minNT(SumF,transpose(x1),var);

Bx&＃61;Funval(FE,var,x2);

if u*Bx

if norm(x2-x1)<&＃61;eps

x&＃61;x2;

break;

else

u&＃61;v*u;

x1&＃61;x2;

end

else

if norm(x2-x1)<&＃61;eps

x&＃61;x2;

break;

else

u&＃61;v*u;

x1&＃61;x2;

end

end

end

minf&＃61;Funval(f,var,x);

format short;

syms x1 x2 r1;

>> p&＃61;taylor(x1&＃43;x2-r1*(1/(x1^2-x2)-1/x1),[x1 x2],[0.001 0.002],&＃39;Order&＃39;,3);

>> x&＃61;[0.1;0.2];k&＃61;0.1;e&＃61;0.01;r&＃61;1;[x,minf]&＃61; minNF(p,x,k,r,e)

x &＃61;

0.0015

minf &＃61;

0.0030

Ans&＃61;

0.0045

室友的&＃xff1a;

2-13&＃xff1a;

内点&＃xff1a;

惩罚函数&＃xff1a;

function anll&＃61;neicheng(p,x,k,r,e)

syms x1 x2 r1;

flag1&＃61;1;

while (flag1)

pd&＃61;subs(p,r1,r);

xold&＃61;x;

flag2&＃61;1;

while (flag2)

dp&＃61;-[diff(pd,x1);diff(pd,x2)];

h&＃61;hessian(pd,[x1,x2]);

h1&＃61;h^-1;

dp_temp&＃61;subs(dp,x1,x(1));

dp_temp&＃61;subs(dp_temp,x2,x(2));

nor&＃61;norm(dp_temp);

if(nor>&＃61;e)

x&＃61;x&＃43;h1*dp_temp;

else

flag2&＃61;0;

end

end

x_temp&＃61;x;

nor2&＃61;norm(x_temp-xold);

if double(nor2)>&＃61;e

r&＃61;k*r;

else

flag1&＃61;0;

end

end

anll&＃61;double(x);

结果&＃xff1a;

clear all; syms x1 x2 r1;

>> p&＃61;taylor(x1&＃43;x2-r1*(1/(x1^2-x2)-1/x1),[x1 x2],[0.001 0.002],&＃39;Order&＃39;,3);

>> x&＃61;[0.1;0.2];k&＃61;0.1;e&＃61;0.01;r&＃61;1;anll&＃61;neicheng(p,x,k,r,e)

anll &＃61;

0.0015

0.0030

Ans&＃61;

0.0045

外点&＃xff1a;

惩罚函数&＃xff1a;

function annn&＃61;waicheng(p,x,k,r,e)

syms x1 x2 r1;

flag1&＃61;1;

while (flag1)

pd&＃61;subs(p,r1,r);

xold&＃61;x;

flag2&＃61;1;

while (flag2)

dp&＃61;-[diff(pd,x1);diff(pd,x2)];

h&＃61;hessian(pd,[x1,x2]);

h1&＃61;h^-1;

dp_temp&＃61;subs(dp,x1,x(1));

dp_temp&＃61;subs(dp_temp,x2,x(2));

nor&＃61;norm(dp_temp);

if(nor>&＃61;e)

x&＃61;x&＃43;h1*dp_temp;

else

flag2&＃61;0;

end

end

x_temp&＃61;x;

nor2&＃61;norm(x_temp-xold);

if double(nor2)>&＃61;e

r&＃61;k*r;

else

flag1&＃61;0;

end

end

annn&＃61;double(x);

结果&＃xff1a;

clear all; syms x1 x2 r1;

p&＃61;taylor(x1&＃43;x2,[x1 x2],[0.001 0.002],&＃39;Order&＃39;,3);

x&＃61;[0.1;0.2];k&＃61;0.1;e&＃61;0.01;r&＃61;1;annn&＃61;waicheng(p,x,k,r,e)

annn &＃61;

0.0015

0.0030

Ans&＃61;

0.0045

2-14

我的(貌似这题抄的他的)&＃xff1a;

2-14

function [x,minf] &＃61; minMixFun(f,g,h,x0,r0,c,var,eps)

gx0 &＃61; Funval(g,var,x0);

if gx0 >&＃61; 0;

else

disp(&＃39;初始点必须满足不等式约束&＃xff01;&＃39;);

x &＃61; NaN;

minf &＃61; NaN;

return;

end

if r0 <&＃61; 0

disp(&＃39;初始障碍因子必须大于0&＃xff01;&＃39;);

x &＃61; NaN;

minf &＃61; NaN;

return;

end

if c >&＃61; 1 || c <0

disp(&＃39;缩小系数必须大于0且小于1&＃xff01;&＃39;);

x &＃61; NaN;

minf &＃61; NaN;

return;

end

if nargin &＃61;&＃61; 7

eps &＃61; 1.0e-6;

end

FE &＃61; 0;

for i&＃61;1:length(g)

FE &＃61; FE &＃43; 1/g(i);

end

FH &＃61; transpose(h)*h;

x1 &＃61; transpose(x0);

x2 &＃61; inf;

while 1

FF &＃61; r0*FE &＃43; FH/sqrt(r0);

SumF &＃61; f &＃43; FF ;

[x2,minf] &＃61; minNT(SumF,transpose(x1),var);

if norm(x2 - x1)<&＃61;eps

x &＃61; x2;

break;

else

r0 &＃61; c*r0;

x1 &＃61; x2;

end

end

minf &＃61; Funval(f,var,x);

Funval.m

function fv &＃61; Funval(f,varvec,varval)

var &＃61; findsym(f);

varc &＃61; findsym(varvec);

s1 &＃61; length(var);

s2 &＃61; length(varc);

m &＃61;floor((s1-1)/3&＃43;1);

varv &＃61; zeros(1,m);

if s1 ~&＃61; s2

for i&＃61;0: ((s1-1)/3)

k &＃61; findstr(varc,var(3*i&＃43;1));

index &＃61; (k-1)/3;

varv(i&＃43;1) &＃61; varval(index&＃43;1);

end

fv &＃61; subs(f,var,varv);

else

fv &＃61; subs(f,varvec,varval);

end

Syms x1 x2;

f&＃61;x1^2-x2^2-3*x2;

g&＃61;1-x1;

h&＃61;x2-2;

[x,minf]&＃61;minMixFun(f,g,h,[2,2],2,0.5,[x1 x2 ],0.001)

x &＃61;

1.0015

minf&＃61;

2.0002

室友的&＃xff1a;

2-14&＃xff1a;

混合&＃xff1a;

惩罚函数&＃xff1a;

function [x,minf] &＃61; MixPunish(f,g,h,x0,r0,c,var,eps)

gx0 &＃61; Funval(g,var,x0);

if gx0 >&＃61; 0;

else

disp(&＃39;初始点必须满足不等式约束&＃xff01;&＃39;);

x &＃61; NaN;

minf &＃61; NaN;

return;

end

if r0 <&＃61; 0

disp(&＃39;初始障碍因子必须大于0&＃xff01;&＃39;);

x &＃61; NaN;

minf &＃61; NaN;

return;

end

if c >&＃61; 1 || c <0

disp(&＃39;缩小系数必须大于0且小于1&＃xff01;&＃39;);

x &＃61; NaN;

minf &＃61; NaN;

return;

end

if nargin &＃61;&＃61; 7

eps &＃61; 1.0e-6;

end

FE &＃61; 0;

for i&＃61;1:length(g)

FE &＃61; FE &＃43; 1/g(i);

end

FH &＃61; transpose(h)*h;

x1 &＃61; transpose(x0);

x2 &＃61; inf;

while 1

FF &＃61; r0*FE &＃43; FH/sqrt(r0);

SumF &＃61; f &＃43; FF ;

[x2,minf] &＃61; minNT(SumF,transpose(x1),var);

if norm(x2 - x1)<&＃61;eps

x &＃61; x2;

break;

else

r0 &＃61; c*r0;

x1 &＃61; x2;

end

end

minf &＃61; Funval(f,var,x);

Funval.m

function fv &＃61; Funval(f,varvec,varval)

var &＃61; findsym(f);

varc &＃61; findsym(varvec);

s1 &＃61; length(var);

s2 &＃61; length(varc);

m &＃61;floor((s1-1)/3&＃43;1);

varv &＃61; zeros(1,m);

if s1 ~&＃61; s2

for i&＃61;0: ((s1-1)/3)

k &＃61; findstr(varc,var(3*i&＃43;1));

index &＃61; (k-1)/3;

varv(i&＃43;1) &＃61; varval(index&＃43;1);

end

fv &＃61; subs(f,var,varv);

else

fv &＃61; subs(f,varvec,varval);

end

Syms x1 x2;

f&＃61;x1^2-x2^2-3*x2;

g&＃61;1-x1;

h&＃61;x2-2;

[x,minf]&＃61;MixPunish(f,g,h,[2,2],2,0.5,[x1 x2 ],0.001)

x &＃61;

1.0015

minf&＃61;

2.0002

结束语

无聊到这地步想必也是没谁了。不过&＃xff0c;刚考完&＃xff0c;我总不能一直玩手机啊。前几天重新看《盘龙》让我在手机上刚了一星期&＃xff0c;不能再这么毫无节制的玩耍了&＃xff0c;但是又不想学习&＃xff0c;所以只好写简书了~~

不过网上毕竟这方面的资源不是很多&＃xff0c;我就算是为后来人做点好事吧&＃xff0c;让你们好找一点~~

个人宣言

知识传递力量&＃xff0c;技术无国界&＃xff0c;文化改变生活&＃xff01;