Set up a gradient descent algorithm for a multivariable economic dispatch problem

Hi,

I'm working on an electrical engineering optimization problem. I tried running the problem following an algorithm which was given to me in class, but I am getting some undefined results. I am attaching my code below:

clc; clear;

syms x y z

vX = [x y z].'; % set up x vector and transposition

Q = [1.562E-3 0 0; 0 1.94E-3 0; 0 0 4.82E-3]; b = [7.92 7.85 7.97]'; % calculate Q and b values according to quadratic problem formula, reproduced as f below

tau = 0.00001; % set up tolerance value

f = vX.'*Q*vX+b.'*vX+949;

Delf = [diff(f,x); diff(f,y); diff(f,z)];

t = (Delf.'*Delf)/(Delf.'*Q*Delf);

vx_old = [300 300 100]'; % set up initial vX value

grad = double(subs(Delf,vX,vx_old));

lenGrad = norm(grad);

Data = [vx_old;lenGrad];

while lenGrad > tau

t = subs(t,vX,vx_old);

vx_new = vx_old-t*grad;

grad = double(subs(Delf,vX,vx_new));

lenGrad = norm(grad);

Data = [Data [vx_new;lenGrad]];

vx_old = vx_new;

end

vx_new = double(vx_new)

fval = double(subs(f,vX,vx_new));

disp(['Minimum value of the function is: ',num2str(fval)]);

disp(['argmin of x value is: ',num2str(vx_new(1))]);

disp(['argmin of y value is: ',num2str(vx_new(2))]);

disp(['argmin of z value is: ',num2str(vx_new(3))]);

For reference, the objective function, which I am trying to minimize, is:

0.001562P1^2 + 0.00194P2^2 + 0.00482P3^2 + 7.92P1 + 7.85P2 + 7.97P3 + 949,

where each P represents the power output of 3 generating thermal units. I understand other methods are better suited to this problem, but I am required to try out unconstrained optimization on this problem. Would you help me determine if my code is wrong at some point or what I can do to improve it?

Answer 1

[Torsten]

There must be something wrong in the computation of t in your code. For t=1, the usual gradient descent works:

clc; clear;
syms x y z 
vX = [x y z].'; % set up x vector and transposition
Q = [1.562E-3 0 0; 0 1.94E-3 0; 0 0 4.82E-3]; b = [7.92 7.85 7.97]'; % calculate Q and b values according to quadratic problem formula, reproduced as f below
tau = 0.00001; % set up tolerance value
f = vX.'*Q*vX+b.'*vX+949;
Delf = [diff(f,x); diff(f,y); diff(f,z)];
sol = solve([diff(f,x)==0, diff(f,y)==0, diff(f,z)==0],[x,y,z]);
double(sol.x)

ans = -2.5352e+03

double(sol.y)

ans = -2.0232e+03

double(sol.z)

ans = -826.7635

t = (Delf.'*Delf)/(Delf.'*Q*Delf);
vx_old = [300 300 100]'; % set up initial vX value
grad = double(subs(Delf,vX,vx_old));
lenGrad = norm(grad);
Data = [vx_old;lenGrad];
iter = 0;
while lenGrad > tau & iter < 10000
    iter = iter + 1;
    %ti = double(subs(t,vX,vx_old));
    %vx_new = vx_old-ti*grad;
    vx_new = vx_old-grad;
    grad = double(subs(Delf,vX,vx_new));
    lenGrad = norm(grad);
    Data = [Data [vx_new;lenGrad]];
    vx_old = vx_new;
end
iter

iter = 4377

vx_new = double(vx_new)

vx_new = 3x1

1.0e+03 * -2.5352 -2.0232 -0.8268

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = double(subs(f,vX,vx_new));
disp(['Minimum value of the function is: ',num2str(fval)]);

Minimum value of the function is: -20326.1329

disp(['argmin of x value is: ',num2str(vx_new(1))]);

argmin of x value is: -2535.2081

disp(['argmin of y value is: ',num2str(vx_new(2))]);

argmin of y value is: -2023.1958

disp(['argmin of z value is: ',num2str(vx_new(3))]);

argmin of z value is: -826.7635

[Mario]

Thank you for your response :) I have some doubts because, in your code, you include two expressions for vx_new on the while loop, and no vx_old update is made. I think I also get different values than what you get on here

[Torsten]

The result agrees with the analytical solution.

I don't understand what you mean with : you include two expressions for vx_new on the while loop, and no vx_old update is made. I just changed

ti = double(subs(t,vX,vx_old));
vx_new = vx_old-ti*grad;

(which would be the correct update of vx_old with your step-size rule) to

vx_new = vx_old-grad;

in your code.

[Mario]

Torsten,

Thank you, and I think I just misunderstood your code at first haha. I did have a few doubts remaining. First, since I am only substracting the vx value by the gradient in your code (within the while loop), does that mean that I do not use the defined step size, t, anymore? Second, since I get a negative value for this problem, which does not make physical sense, would you have any advice on how to improve the code to better reflect thermal generating units operating conditions?

[Torsten]

First, since I am only substracting the vx value by the gradient in your code (within the while loop), does that mean that I do not use the defined step size, t, anymore?

The step size is t=1 in each step.

Second, since I get a negative value for this problem, which does not make physical sense, would you have any advice on how to improve the code to better reflect thermal generating units operating conditions?

Your f gives a negative value for its minimum - thus you have to modify f. Since I don't understand the background of your problem, I cannot give further advice.

By the way: Your t is not correct; it must read

t = (Delf.'*Delf)/(Delf.'*(2*Q)*Delf);

Replacing it in your code, convergence is really fast:

clc; clear;
syms x y z 
vX = [x y z].'; % set up x vector and transposition
Q = [1.562E-3 0 0; 0 1.94E-3 0; 0 0 4.82E-3]; b = [7.92 7.85 7.97]'; % calculate Q and b values according to quadratic problem formula, reproduced as f below
tau = 0.00001; % set up tolerance value
f = vX.'*Q*vX+b.'*vX+949;
Delf = [diff(f,x); diff(f,y); diff(f,z)];
sol = solve([diff(f,x)==0, diff(f,y)==0, diff(f,z)==0],[x,y,z]);
double(sol.x)

ans = -2.5352e+03

double(sol.y)

ans = -2.0232e+03

double(sol.z)

ans = -826.7635

t = (Delf.'*Delf)/(Delf.'*(2*Q)*Delf);
vx_old = [300 300 100]'; % set up initial vX value
grad = double(subs(Delf,vX,vx_old));
lenGrad = norm(grad);
Data = [vx_old;lenGrad];
iter = 0;
while lenGrad > tau & iter < 10000
    iter = iter + 1;
    ti = double(subs(t,vX,vx_old));
    vx_new = vx_old-ti*grad;
    grad = double(subs(Delf,vX,vx_new));
    lenGrad = norm(grad);
    Data = [Data [vx_new;lenGrad]];
    vx_old = vx_new;
end
iter

iter = 20

vx_new = double(vx_new)

vx_new = 3x1

1.0e+03 * -2.5352 -2.0232 -0.8268

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = double(subs(f,vX,vx_new));
disp(['Minimum value of the function is: ',num2str(fval)]);

Minimum value of the function is: -20326.1329

disp(['argmin of x value is: ',num2str(vx_new(1))]);

argmin of x value is: -2535.2091

disp(['argmin of y value is: ',num2str(vx_new(2))]);

argmin of y value is: -2023.1959

disp(['argmin of z value is: ',num2str(vx_new(3))]);

argmin of z value is: -826.763