I am trying to code a steepest descent problem. I have the main file:
I have edited my file and combined the two functions into one .
The objective function is "c". The variables in this problem are q, t and alpha. q and t appear in "c", while alpha appears in the re-expressed anonymous function based off "c" in the later part of the code.
I'm trying to minimize c using the steepest descent problem. First I find the gradient (partial differential equations of c with respect to q and t), then i use subs and double to generate new values of q and t with each iteration in the while loop. I also use fminbnd as a line search to find the value of the step size, "alpha" that minimizes the objective function during each iteration.
syms q t alphaa = 0.20; %annual fixed charges, fraction
c_c = 12.50; %crude oil price, $/kL
c_i = 0.50; %insurance cost, $/kL
c_x = 0.90; %customs cost, $/kL
i = 0.10; %interest rate
n = 2; %number of ports
p = 7000; %land price, $/m2
c = c_c + c_i + c_x + ... + (2.09e4*t.^-0.3017) / 360 ... + (1.064e6*a*t.^0.4925) / (52.47*q*360) ... + (4.242e4*a*t^0.7952) / (52.47*q.*360) ... + (1.813*i*p*(n*t+1.2*q).^0.861) / (52.47*q*360) ... + 4.25e3*a*(n*t+1.2*q) / (52.47*q*360) ... + (5.024e3*q.^-0.1899)/360 + (0.1049*q.^0.671)/360;k=0; % iteration counter start for outer loop
eps=1; %tolerance for outer loop
x=[150000;450000]; %initial starting point
dq=diff(c,q); %equation of partial differentiation of c wrt q
dt=diff(c,t); %equation of partial differentiation of c wrt t
alpha = 100000; %initial step size for outer loop
while eps>1e-10 & k<10000 %tolerance and number of iterations
eps=double(abs(subs(dq, [q,t], [x(1),x(2)]))+abs(subs(dt, [q,t], [x(1),x(2)]))); % I used subs & double to make the partial differentiation terms dq and dt
% scalar to calculate eps.
a1 = x(1) - alpha*double(subs(dq,[q,t],[x(1),x(2)]));a2 = x(2) - alpha*double(subs(dt,[q,t],[x(1),x(2)]));[alpha1,fval] = fminbnd(@(alpha) (c_c + c_i + c_x + ... + (2.09e4*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)]).^-0.3017) / 360 ... + (1.064e6*a*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)]).^0.4925) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])*360) ... + (4.242e4*a*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)])^0.7952) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)]).*360) ... + (1.813*i*p*(n*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)])+1.2*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])).^0.861) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])*360) ... + 4.25e3*a*(n*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)])+1.2*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])*360) ... + (5.024e3*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)]).^-0.1899)/360 + (0.1049*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)]).^0.671)/360))))))))))))))))))))))))))),10,120000)alpha = alpha1;x(1)=a1;x(2)=a2;xk=k+1;endThe problem is that I can get the code working, but the values of q, t and alpha are changing extremely slowly in each iteration.
I'm also having problems with "alpha" – in theory, "alpha" is supposed to change with each iteration, however my alpha value always takes on the maximum value (upper bound) in fminbnd and remains stationary.
If I remove this portion of the code:
[alpha1,fval] = fminbnd(@(alpha) (c_c + c_i + c_x + ... + (2.09e4*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)]).^-0.3017) / 360 ... + (1.064e6*a*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)]).^0.4925) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])*360) ... + (4.242e4*a*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)])^0.7952) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)]).*360) ... + (1.813*i*p*(n*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)])+1.2*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])).^0.861) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])*360) ... + 4.25e3*a*(n*(x(2)-alpha*double(subs(dt,[q,t],[x(1),x(2)])+1.2*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])) / (52.47*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)])*360) ... + (5.024e3*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)]).^-0.1899)/360 + (0.1049*(x(1)-alpha*double(subs(dq,[q,t],[x(1),x(2)]).^0.671)/360))))))))))))))))))))))))))),10,120000)alpha = alpha1;The convergence rate increases greatly. Any idea as to what may be causing this issue?
Best Answer