- Utilizing The Code Provided In Matlab Code Function X0 Lambda Y Main X0 Lambda Y Mu If Nargin 4 X0 40 40 40 4 1 (71.71 KiB) Viewed 28 times
function [x0,lambda,y]=main(x0,lambda,y,mu)
if nargin<4
x0=[40;40;40;-40];
lambda=[0;0];
y=[0;0];
mu=200;
end
format compact
Ae=[1,1,1,1;-1,1,1,1];
be=[2;3];
Ai=[1,-1,1,1;1,1,-1,1];
bi=[1;4];
f=@(x) norm(x)^2
%x0=[40;40;40;-40]
h=@(x) Ae*x-be;
g=@(x) Ai*x-bi;
%mu=200;
s=-g(x0)
%lambda=[0;0];
%y=[0;0];
for k=1:400
[A,B]=fullsystem(f,g,h,x0,s,lambda,y,mu);
ds=-A\B;
x0=x0+ds(1:4);
s=s+inv(diag(s))*ds(5:6);
y=y+ds(7:8);
lambda=lambda+ds(9:10);
R=norm(ds([1:4]));
if R<1e-5
break
end
end
x0
Ae*x0-be
s
end
function [A,B]=fullsystem(f,g,h,x,s,lambda,y,mu)
H=Hessian(f,g,h,x,lambda,y);
Jg=Jacobiang(g,x);
Jh=Jacobianh(h,x);
nx=length(x);
ng=length(s);
nh=length(h(x));
A=[H,zeros(nx,ng),Jh' ,Jg' ;zeros(ng,nx),
diag(s)*diag(lambda),zeros(ng,nh),diag(s);Jh,zeros(nh,2*ng+nh);Jg,diag(s),zeros(ng,nh+ng)];
B=[gradientf(f,x)+Jh'*y+Jg'*lambda;diag(s)*lambda-mu*ones(ng,1);h(x);g(x)+s];
end
function gf=gradientf(f,x)
dx=1e-5;
gf=zeros(length(x),1);
for j=1:length(x)
dxj=0*x;
dxj(j)=dx;
gf(j)=(f(x+dxj)-f(x))/dx;
end
end
function H=Hessian(f,g,h,x,lambda,y)
dx=1e-5;
H=zeros(length(x),length(x));
for j=1:length(x)
dxj=0*x;
dxj(j)=dx;
dfj= @(x) (f(x+dxj)-f(x))/dx;
for k=j:length(x)
dxk=0*x;
dxk(k)=dx;
H(j,k)=(dfj(x+dxk)-dfj(x))/dx;
H(k,j)=H(j,k);
end
end
G=zeros(length(x),length(x));
for j=1:length(x)
dxj=0*x;
dxj(j)=dx;
dgj=@(x) (g(x+dxj)-g(x))/dx;
for k=j:length(x)
dxk=0*x;
dxk(k)=dx;
ddgjk=(dgj(x+dxk)-dgj(x))/dx;
G(j,k)=sum(ddgjk.*lambda);
G(k,j)=G(j,k);
end
end
Y=zeros(length(x),length(x));
for j=1:length(x)
dxj=0*x;
dxj(j)=dx;
dhj=@(x) (h(x+dxj)-h(x))/dx;
for k=j:length(x)
dxk=0*x;
dxk(k)=dx;
ddhjk=(dhj(x+dxk)-dhj(x))/dx;
Y(j,k)=sum(ddhjk.*y);
Y(k,j)=Y(j,k);
end
end
H=H+G+Y;
end
function Jg=Jacobiang(g,x)
dx=1e-5;
Jg=zeros(length(g(x)),length(x));
for j=1:length(x)
dxj=0*x;
dxj(j)=dx;
Jg(:,j)=(g(x+dxj)-g(x))/dx;
end
end
function Jh=Jacobianh(h,x)
dx=1e-5;
Jh=zeros(length(h(x)),length(x));
for j=1:length(x)
dxj=0*x;
dxj(j)=dx;
Jh(:,j)=(h(x+dxj)-h(x))/dx;
end
end
As demonstrated in class (specifically see lectures 4/15, 4/18, and 4/20) implement an interior point method similar to the interior point algorithm in fmincon except only with a Newton step (no-CG step is required). Use it to solve the following constrained minimization problems: 1) min x} + x2 + x3 + x2 + xz subject to X1 – x2 + x3 = 3 and x3 + x4 > 3 and 25 < -2 2) min x1 + x2 + xj subject to sin(x1 + x2 + x3) = V(2)/2 and xz + xz> 1 + Include all your code and a table showing approximately 20-30 evenly spaced iterations and the values of x, s, lambda, y and mu. (So if your algorithm ran in 250 iterations, include every 10th iteration in the table, if it runs in 5000, include every 200th iteration in the table)