Power loss minimisation using fmincon

I need to minimise my power losses using fmincon with constraints.

the Values that i am using are:

Vm=230; Rmk=0.0246; Rkn=0.048; Gmk = 33.774; Bmk = -15.240; Gkn=19.745; Bkn=-4.628;

My objective function is:

P_loss = @(V) Rmk*abs(Ymk*(Vm-Vk))^2+Rkn*abs(Ykn*(Vk-Vn))^2;

My power flow equations that i have are:

function [c,ceq] = power_flow_equations(V)
c = [];
Pn_EV = 4000; Vmx = 230; Gmk = 33.774; Bmk = -15.240; Gkn=19.745; Bkn=-4.628; Pk=59280; Pn_others=53040; Qn_others = 32870; Qk=19480;
ceq(1) = V(1)^2*(Gmk+Gkn)+V(2)^2*(Gmk+Gkn)-V(1)*Vmx*Gmk-V(2)*Vmx*Bmk-V(1)*V(3)*Gkn+V(1)*V(4)*Bkn-V(2)*V(3)*Bkn-V(2)*V(4)*Gkn+Pk;
ceq(2) = V(1)^2*(Bmk+Bkn)+V(2)^2*(Bmk+Bkn)-V(1)*Vmx*Bmk+V(2)*Vmx*Gmk-V(1)*V(3)*Bkn-V(1)*V(4)*Gkn+V(2)*V(3)*Gkn-V(2)*V(4)*Bkn-Qk;
ceq(3) = Gkn*V(3)^2+Gkn*V(4)^2-Gkn*V(3)*V(1)-Gkn*V(4)*V(2)+Bkn*V(3)*V(2)-Bkn*V(4)*V(1)+Pn_others+Pn_EV;
ceq(4) = Bkn*V(3)^2+Bkn*V(4)^2-Gkn*V(3)*V(2)+Gkn*V(4)*V(1)-Bkn*V(3)*V(1)-Bkn*V(4)*V(2)-Qn_others;
end

I want to add the voltage constraint which is:

216.2<=sqrt(Vx^2+Vy^2)<=253

and add some power constraints like

-4kW<=Pn_EV<=4kW

This is what is what i have so far

clc
clear
Vm=230; Rmk=0.0246; Rkn=0.048; Gmk = 33.774; Bmk = -15.240; Gkn=19.745; Bkn=-4.628;
V=[150*ones(4,1)];
Vk=V(1)+j*V(2);
Vn=V(2)+j*V(4);
Ymk=Gmk+j*Bmk;%Representing the Admitance from m to k
Ykn=Gkn+j*Bkn;%Representing the Admitance from k to n   
format shortG makes the numbers nicer to look at for me
e in Matlab is a shorthand for *10^
format longG
P_loss = @(V) Rmk*abs(Ymk*(Vm-Vk))^2+Rkn*abs(Ykn*(Vk-Vn))^2;
A=[];
b=[];
Aeq = [2*V(1)*(Gkn + Gmk) + Bkn*V(4) - Gkn*V(3) - Gmk*Vm, 2*V(2)*(Gkn + Gmk) - Bkn*V(3) - Bmk*Vm - Gkn*V(4), - Bkn*V(2) - Gkn*V(1), Bkn*V(1) - Gkn*V(2)
       2*V(1)*(Bkn + Bmk) - Bkn*V(3) - Bmk*Vm - Gkn*V(4), 2*V(2)*(Bkn + Bmk) - Bkn*V(4) + Gkn*V(3) + Gmk*Vm, Gkn*V(2) - Bkn*V(1), - Bkn*V(2) - Gkn*V(1)
       - Bkn*V(4) - Gkn*V(3), Bkn*V(3) - Gkn*V(4), Bkn*V(2) - Gkn*V(1) + 2*Gkn*V(3), 2*Gkn*V(4) - Gkn*V(2) - Bkn*V(1)
       Gkn*V(4) - Bkn*V(3), - Bkn*V(4) - Gkn*V(3), 2*Bkn*V(3) - Bkn*V(1) - Gkn*V(2), 2*Bkn*V(4) - Bkn*V(2) + Gkn*V(1)];
beq=zeros(4,1);
Aeq=[];
beq=[];
nonlcon= @power_flow_equations;

lb = [sqrt(216^2-V(2)^2);sqrt(216^2-V(1)^2);sqrt(216^2-V(4)^2);sqrt(216^2-V(3)^2)];
ub = [sqrt(253^2-V(2)^2);sqrt(253^2-V(1)^2);sqrt(253^2-V(4)^2);sqrt(253^2-V(3)^2)];
lb = [];
ub = [];
options.Display = 'iter';
[x,fval,exitflag]=fmincon(P_loss,V,A,b,Aeq,beq,lb,ub,nonlcon,options)

Vk=V(1)+j*V(2);
Vn=V(2)+j*V(4);
exitflag
formatSpec = 'Power loss=%.3f watts\n';
fprintf(formatSpec,P_loss)
disp('Solution')
formatSpec = 'Vkx=%.3f Vky=%.3f, Vnx=%.3f, Vny=%.3f Vk=%.3f Vn=%.3f\n';
fprintf(formatSpec, x(1),x(2),x(3),x(4),abs(Vk),abs(Vn))

I do not know how i can implement my power contraints into this or set the upper and lower voltage bounds. Any help will be appreciated. The unknown variables are V(1) V(2) V(3) and V(4) and the power of the EV Pn_EV

[John D'Errico]

When you delete your question, you insult everyone who bothered to waste their time in answering your question. This is rude to them. It is harmful to the site, because now your question is completely useless to anyone else who might have been able to learn from their answers.

[Rik]

Are you trying to see if you can be more stubborn than me? I will keep restoring your question until someone from Mathworks tells me otherwise. You could also try explaining why you want to delete this question in the first place.

[William Rose]

@John D'Errico and @Rik,

Thank you for maintaining the integrity of the online dicussion. Your efforts benefit the users of this site.

Answer 1

[William Rose]

Unlike some people who post on this site, you have done a good job of explaining your problem and attempting at a solution and posting your attempt at a solution.

I have some suggestions and questions.

1. I think you need to adjust P_loss() for this to work correctly. Specifically, your code includes

Vk=V(1)+j*V(2);
Vn=V(2)+j*V(4);

and later it has

P_loss = @(V) Rmk*abs(Ymk*(Vm-Vk))^2+Rkn*abs(Ykn*(Vk-Vn))^2;

When fmincon() tries to minimize P_loss() by adjusting the four elements of V, it will use the old (initial) values of Vk and Vn in its calculations - but it should use the current trial values of V(:). So you should rewrite P_loss using V(1), V(2), V(3), V(4) intead of Vk and Vn. Or make P_loss a separate function which computes Vk and Vn internally before it compute P_loss.

2. You say tha you have a voltage constraint:

216.2<=sqrt(Vx^2+Vy^2)<=253

What are Vx and Vy?

3. You say you have power constraint:

-4kW<=Pn_EV<=4kW

But your function power_flow_equations() defines

    Pn_EV=4000;

I do not undertand how Pn_EV is a constraint, since it is defined as a constant..

4. Your first two equations for lb and ub are moot, because the next two equations undo the first two. Which I think you know.

[William Rose]

@Lauris Cuzans

OK. What is the formula for Pn_EV? I assume it is a function of the voltage, and maybe of the constants.

[William Rose]

Your changes to power_flow() and to P_loss() look good. P_loss(), as you have defined it, does not depend at all on v(5). Thefore I predict fmincon() will search unsuccessfully for a minimum and will quit when the maximum number of iterations is reached. It will probably issue a warning when it finishes, to say that it was unable to find a minimum.

The constraint -4000<=Pn_EV<=4000 seems to be a formulation of equations 12 and 13 of the paper you sent, where Pmax=4000, and

. ,

where your Pn_EV is on the left hand side above, and equivalent terms from the paper are on the right. But I don't know how the right hand side terms relate to the variable definitions in your code.

Even if you do develop an equation that defines Pn_EV in terms of V and the constants, fmincon() will not be able to find a minimum, as long as the P_loss() does not depend on V(5).

Answer 2

[William Rose]

I would make power_flow_equations() a standalone function and pass the reference to it to fmincon(). Since you do not have specific direct contraints on V(:), at least not that your have mentioned, I would replace lb and ub with [] in the call to fmincon().

[x,fval,exitflag]=fmincon(P_loss,V,A,b,Aeq,beq,[],[],@power_flow_equations,options);

To impose a power constraint (once you define Vx and Vy, as I mentioned in my previous post), you can add two inequality constraints to power_flow_equations(), which will be elements 1 and 2 of vector c. Therefore the edited power_flow_equations() will be something like this:

function [c,ceq] = power_flow_equations(V)
%add code to define Vx, Vy in terms of V()
c(1) = sqrt(Vx^2+Vy^2)-253;     %fmincon will enforce c(1)<=0 i.e. sqrt(Vx^2+Vy^2)<=253
c(2) = 216.2-sqrt(Vx^2+Vy^2);   %fmincon will enforce c(2)<=0 i.e. 216.2<=sqrt(Vx^2+Vy^2)
%include the code for the equality constraints, as before
Pn_EV = 4000; Vmx = 230; Gmk = 33.774; Bmk = -15.240; Gkn=19.745; Bkn=-4.628; Pk=59280; Pn_others=53040; Qn_others = 32870; Qk=19480;
ceq(1) = V(1)^2*(Gmk+Gkn)+V(2)^2*(Gmk+Gkn)-V(1)*Vmx*Gmk-V(2)*Vmx*Bmk-V(1)*V(3)*Gkn+V(1)*V(4)*Bkn-V(2)*V(3)*Bkn-V(2)*V(4)*Gkn+Pk;
ceq(2) = V(1)^2*(Bmk+Bkn)+V(2)^2*(Bmk+Bkn)-V(1)*Vmx*Bmk+V(2)*Vmx*Gmk-V(1)*V(3)*Bkn-V(1)*V(4)*Gkn+V(2)*V(3)*Gkn-V(2)*V(4)*Bkn-Qk;
ceq(3) = Gkn*V(3)^2+Gkn*V(4)^2-Gkn*V(3)*V(1)-Gkn*V(4)*V(2)+Bkn*V(3)*V(2)-Bkn*V(4)*V(1)+Pn_others+Pn_EV;
ceq(4) = Bkn*V(3)^2+Bkn*V(4)^2-Gkn*V(3)*V(2)+Gkn*V(4)*V(1)-Bkn*V(3)*V(1)-Bkn*V(4)*V(2)-Qn_others;
end

Try.

Answer 3

[William Rose]

You said that one of the constraints you would like to impose is

-4kW<=Pn_EV<=4kW

but then you defined Pn_EV as a constant, which does not make sense. Maybe that was a temporary measure to get your code working. If Pn_EV is actually a function of V(), then you can impose the desired constraint by adding c(3) and c(4) to the code for power_flow_equations(), as shown below. I will assume that Pn_EV=(V(1)^2-V(2)^2)/Gmk. You should replace that with the correct equation for Pn_EV.

function [c,ceq] = power_flow_equations(V)
%Define the constants first. Delete Pn_EV from the list of constants.
Vmx = 230; Gmk = 33.774; Bmk = -15.240; Gkn=19.745; Bkn=-4.628; 
Pk=59280; Pn_others=53040; Qn_others = 32870; Qk=19480;
%Add code to define Vx, Vy in terms of V().
%Define Pn_EV in terms of V() and the constants.
Pn_EV=(V(1)^2-V(2)^2)/Gmk;      %adjust as needed
c(1) = sqrt(Vx^2+Vy^2)-253;     %fmincon will enforce c(1)<=0 i.e. sqrt(Vx^2+Vy^2)<=253
c(2) = 216.2-sqrt(Vx^2+Vy^2);   %fmincon will enforce c(2)<=0 i.e. 216.2<=sqrt(Vx^2+Vy^2)
c(3) = Pn_EV-4000;              %fmincon will enforce c(3)<=0 i.e. Pn_EV<=4000
c(4) = -4000-Pn_EV;             %fmincon will enforce c(4)<=0 i.e. -4000<=Pn_EV
%and so on for ceq()
end

Try.