Display of result in tabular form
1 view (last 30 days)
Show older comments
Good morning friends, Am trying to solve an unconstrained optimization problem using (trust region) dogleg method. Below is a program i wrote for the problem, please i need help on how to make the results be in a tabular form under their respective headings and also plotting the trust region radius at each point. My program; clear all syms x1 x2 X = [x1 ; x2]; f = (x2 - 0.129*x1^2 + 1.6*x1 - 6)^2 + 6.07*cos(x1) + 10; w = jacobian(f,X); g1 = w'; h1 = jacobian(jacobian(f,X)); k = 0; x1 = 6; x2 = 14; I = 2.0; M = 5.0; B = 0.25; C = 2.0; lambda = 0.01; g = eval(g1); h = eval(h1); b = norm(g); while k<20 [V,D] = eig(h); [m,n] = size(h); d = diag(D); dmin = min(d); for i=1:m subh=h(1:i,1:i); if (det(subh)>0) H = h; else H = h + eye(size(h))*(lambda-dmin); end end Y = inv(H); r = (g'*g)/(g'*H*g); pU = -r*g; pB = -Y*g; if 0<=r<=1 P = r*pU; elseif 1<=r<=2 P = pU + (r-1)*(pB-pU); end PP = norm(P); b = norm(g); %where b is tol
f = @(x1,x2) (x2-0.129*x1^2+1.6*x1-6)^2+6.07*cos(x1)+10;
X = [x1;x2];
f0 = f(x1,x2);
W = g'*P;
Z = 0.5*(P'*H*P);
y = f0+W+Z;
v = f0;
A = X+P;
x1 = A(1,1);
x2 = A(2,1);
q = f(x1,x2);
Q = 0.2; R = 0.25; S = 0.75;
e = (f0 - q)/(v - y);
if e<R
Inew = B*I;
elseif e>S && PP==I
Inew = min(C*I,M);
else
Inew = I;
end
if e>Q
Xnew = A;
else
Xnew = X;
end
disp( [k' f0' x1' x2' e' I' P' b' PP'] )
Xnew = subs(Xnew);
k = k+1;
X = Xnew;
x1 = X(1,1);
x2 = X(2,1);
gnew = eval(g1);
b = norm(gnew);
hnew = eval(h1);
I = Inew; g = gnew;
[Vnew,Dnew] = eig(hnew);
[mnew,nnew] = size(hnew);
dnew = diag(Dnew);
dnewmin = min(dnew);
for i=1:mnew
subhnew=hnew(1:i,1:i);
if (det(subhnew)>0)
Hnew = hnew;
else
Hnew = hnew + eye(size(hnew))*(lambda-dnewmin);
end
end
Ynew = inv(Hnew);
rnew = (gnew'*gnew)/(gnew'*Hnew*gnew);
pUnew = -rnew*gnew;
pBnew = -Ynew*gnew;
if 0<=rnew<=1
Pnew = rnew*pUnew;
elseif 1<=rnew<=2
Pnew = pUnew + (rnew-1)*(pBnew-pUnew);
end
PPnew = norm(Pnew);
bnew = norm(gnew);
V = Vnew; D = Dnew; h = hnew; m = mnew; n = nnew; d = dnew; H = Hnew;
Y = Ynew; r = rnew; pU = pUnew; pB = pBnew; PP = PPnew; P = Pnew;
end
Thanks you as i await your kind help, inputs and corrections.
0 Comments
Answers (1)
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!