Dear Ray et al.
The last script I have posted here had some problem for calculating qji ,
qij. I have modified the code and got almost similar results compare to
MATPOWER. However there are still some small differences in results. Where
is this small difference in results originated?
Here I post the modified version plus results:
Thank you for your consideration .
% calculating line loading via Matpower:
define_constants; mpc = loadcase('case4gs');
RES=runpf(mpc);
flow_from = sqrt(RES.branch(:, PF).^2 + RES.branch(:, QF).^2);
flow_to = sqrt(RES.branch(:, PT).^2 + RES.branch(:, QT).^2);
flow = (flow_from + flow_to)/2;
MATPOWER_loading = flow ./ RES.branch(:, RATE_A);
%%% Here is my script . %%%
% obtained from Matpower case4gs.m
Base=100;from=[1 1 2 3]'; to=[2 3 4 4]';
r=[ 0.01008 0.00744 0.00744 0.01272]';
x=[ 0.0504 0.0372 0.0372 0.0636]';
B=[ 0.1025 0.0775 0.0775 0.1275]';
Rating=[250 250 250 250]';
% Load flow results
v_m=[1 0.982 0.969 1.02]';
v_deg=[0 -0.976 -1.872 1.523]';
v_a = v_deg.*(pi/180); % convert degree to rad
Z_m=sqrt(r.^2+x.^2); Z_phi=acos(r./Z_m);
pij=Base*(v_m(from)./Z_m).*(v_m(from).*cos(Z_phi)-v_m(to).*cos(Z_phi-v_a(to)+v_a(from)));
qij=Base*v_m(from).*((-0.5*B.*v_m(from))+(Z_m.^-1).*(v_m(from).*sin(Z_phi)-v_m(to).*sin(Z_phi-v_a(to)+v_a(from))));
Sij=sqrt(pij.^2+qij.^2);
pji=Base*(v_m(to)./Z_m).*(v_m(to).*cos(Z_phi)-v_m(from).*cos(Z_phi-v_a(from)+v_a(to)));
qji=Base*v_m(to).*((-0.5*B.*v_m(to))+(Z_m.^-1).*(v_m(to).*sin(Z_phi)-v_m(from).*sin(Z_phi-v_a(from)+v_a(to))));
Sji=sqrt(pji.^2+qji.^2);
S=Base*(Sji+Sij)./2;
My_loading=S./(Base*Rating);
table(from,to,MATPOWER_loading,My_loading)
table(pij,RES.branch(:, PF),pji,RES.branch(:, PT), qij,RES.branch(:,
QF),qji,RES.branch(:, QT))
Results:
from to MATPOWER_loading My_loading
____ __ ________________ __________
1 2 0.18842 0.19067
1 3 0.46338 0.46339
2 4 0.60769 0.61037
3 4 0.47707 0.47707
On Thu, Apr 23, 2015 at 2:40 AM, Electric <[email protected]>
wrote:
> Dear Ray et al
> Thank you for your earlier answer.
> Returning to my original question, Implementing line loading constraints
> in for own method, I have tried to obtain the line loading equations and
> compare them with MATPOWER results. To obtain the line loading, I first run
> a power flow calculation, using MATPOWER to get the voltage magnitude and
> angle at each buses.
> Here I post my script. If you have some time, please help me to find the
> reason of this difference:
> % calculating line loading via Matpower:
> define_constants;
> mpc = loadcase('case4gs');
> RES=runpf(mpc);
> flow_from = sqrt(RES.branch(:, PF).^2 + RES.branch(:, QF).^2);
> flow_to = sqrt(RES.branch(:, PT).^2 + RES.branch(:, QT).^2);
> flow = (flow_from + flow_to)/2;
> MATPOWER_loading = flow ./ RES.branch(:, RATE_A);
> % Here is my script .
> % obtained from Matpower case4gs.m
> from=[ 1 1 2 3]';
> to= [ 2 3 4 4]';
> r=[ 0.01008 0.00744 0.00744 0.01272]';
> x=[ 0.0504 0.0372 0.0372 0.0636]';
> B=[ 0.1025 0.0775 0.0775 0.1275]';
> Rating=[ 250 250 250 250]';
> % Load flow results
> v_m=[ 1 0.982 0.969 1.02]';
> v_deg=[ 0 -0.976 -1.872 1.523]';
> v_a = v_deg.*(pi/180); % convert degree to rad
> Z_m=sqrt(r.^2+x.^2);
> Z_phi=acos(r./Z_m);
>
> pij=(v_m(from)./Z_m).*(v_m(from).*cos(Z_phi)-v_m(to).*cos(Z_phi-v_a(to)+v_a(from)));
>
> qij=(v_m(from).*((-0.5*v_m(from).*B)+(sin(Z_phi).*v_m(from)./Z_m)-v_m(to).*sin(Z_phi-v_a(to)+v_a(from))));
> Sij=sqrt(pij.^2+qij.^2);
>
>
> pji=(v_m(to)./Z_m).*(v_m(to).*cos(Z_phi)-v_m(from).*cos(Z_phi-v_a(from)+v_a(to)));
>
> qji=(v_m(to).*((-0.5*v_m(to).*B)+(sin(Z_phi).*v_m(to)./Z_m)-v_m(from).*sin(Z_phi-v_a(from)+v_a(to))));
> Sji=sqrt(pji.^2+qji.^2);
> S=(Sji+Sij)./2;
> My_loading=S./Rating;
> table(from,to,MATPOWER_loading,My_loading)
>
> Output :
> from to MATPOWER_loading My_loading
> ____ __ ________________ __________
>
> 1 2 0.18842 0.070916
> 1 3 0.46338 0.096366
> 2 4 0.60769 0.099696
> 3 4 0.47707 0.055883
>
> Bests.
>
> On Fri, Apr 17, 2015 at 6:05 PM, Ray Zimmerman <[email protected]> wrote:
>
>> Yes, it is possible, as long as your costs can be put in the form of
>> (6.27)-(6-31) in the User’s Manual
>> <http://www.pserc.cornell.edu/matpower/docs/MATPOWER-manual-5.1.pdf>.
>> And you can simply set the appropriate parameters in the gencost matrix
>> to zero, rather than forming user-defined cost terms to cancel out their
>> contribution to the objective function.
>>
>> It looks like your costs are simple linear functions of up and down
>> generation shift variables. The standard optimization variables contain
>> generation, but not separate up and down shifts in generation from a
>> reference. You can define those by introducing a set of user-defined linear
>> constraints, in the form of (6.25) …
>>
>> Pg – Pup + Pdn = Pref
>>
>> … where Pref is a constant reference dispatch, Pg is part of the existing
>> optimization variable *x*, and Pup and Pdn are new vectors of variables
>> (blocks in *z*) defined by additional columns in the *A* matrix. You
>> also need to set z_min to zero to enforce ...
>>
>> Pup >= 0
>> Pdn >= 0
>>
>>
>> Then, I believe, your costs are purely linear functions of these new
>> user-defined optimization variables, so you can use the simplest form of
>> the user-defined cost, where parameters d_i and m_i are all 1, k_i and
>> r_hat are all zero (so *w* = *r*), the *H* matrix is zero and the *C*
>> vector is all ones. In that case, (6.27) and following simply reduces to …
>>
>> f_u(x, z) = C' * N * [x; z]
>>
>> All you need to do then is create N in two blocks, stacked vertically,
>> one for your f1 and one for your f2.
>>
>>
>> --
>> Ray Zimmerman
>> Senior Research Associate
>> B30 Warren Hall, Cornell University, Ithaca, NY 14853 USA
>> phone: (607) 255-9645
>>
>> On Apr 16, 2015, at 1:04 PM, Electric <[email protected]>
>> wrote:
>>
>> Dear Prof Ray
>>
>> I studied extended OPF function to adopt my problem to be run with
>> MATPOWER package as you recommended already, so I could avoid implementing
>> standard OPF constraints like line loadings and voltage range. As I am
>> running a multi-objective congestion management problem (equivalently,
>> maybe called multi-objective OPF), is it possible to include other
>> objective functions separately in the OPF cost function?
>> for example :
>>
>> f1= I have my own cost function, so I can minus the existed cost function
>> to cancel it and add my own cost functions that are still directly depended
>> on the optimization vector. f1=sum(X[i].Bup[i]+Y[i].Bdown[i]), where X and
>> Y are up and down generation shift of unit i and Bup[i] and Bdown[i] are
>> their corresponding bid.
>>
>> f2= security function.
>> f2 is also directly depended on the optimization vector.
>> f2=sum(X[i].S[i]), where x[i] is up generation shift of unit i and S[i] is
>> its corresponding sensitivity to the stability margin.
>>
>> Your help is much appreciated.
>>
>>
>> On Fri, Apr 10, 2015 at 9:11 PM, Electric <
>> [email protected]> wrote:
>>
>>> Dear Prof Ray and Shiri
>>> Thank you so much for your help. I am going to study opf_consfcn()
>>> <http://www.pserc.cornell.edu//matpower/docs/ref/matpower5.1/opf_consfcn.html>
>>> and find out how to implement the constraints I want.
>>> Kind regards.
>>>
>>> On Fri, Apr 10, 2015 at 5:43 PM, Ray Zimmerman <[email protected]> wrote:
>>>
>>>> Specifically, you might look at how the power flow and line limit
>>>> constraints are implemented in opf_consfcn()
>>>> <http://www.pserc.cornell.edu//matpower/docs/ref/matpower5.1/opf_consfcn.html>.
>>>> However, for a problem that includes the standard OPF constraints I still
>>>> think it would be simpler and cleaner to set up your problem to use the
>>>> extended OPF formulation in MATPOWER, unless of course you have some
>>>> requirement the precludes that option.
>>>>
>>>> Ray
>>>>
>>>>
>>>> On Apr 9, 2015, at 10:40 PM, Abhyankar, Shrirang G. <
>>>> [email protected]> wrote:
>>>>
>>>> http://www.pserc.cornell.edu//matpower/docs/ref/
>>>>
>>>> On Apr 9, 2015, at 4:00 PM, "Electric" <[email protected]>
>>>> wrote:
>>>>
>>>> Dear Shri
>>>> Thanks for your answer. But I don't want to use runopf(). I want to
>>>> write these constraint manually with Ybus and equations I've written. Is
>>>> that possible? How should I start?
>>>>
>>>> Thank you.
>>>>
>>>> On Fri, Apr 10, 2015 at 1:24 AM, Abhyankar, Shrirang G. <
>>>> [email protected]> wrote:
>>>>
>>>>> I don't know what your algorithm does but I think you can do the
>>>>> same, determine optimal generation and load, using MATPOWER's optimal
>>>>> power
>>>>> flow (runopf). runopf supports all the constraints that you've listed.
>>>>>
>>>>>
>>>>> Shri
>>>>>
>>>>> ------------------------------
>>>>>
>>>>> *From:* [email protected] [
>>>>> [email protected]] on behalf of Electric [
>>>>> [email protected]]
>>>>> *Sent:* Thursday, April 09, 2015 3:18 PM
>>>>> *To:* MATPOWER discussion forum
>>>>> *Subject:* How to verify the feasibility of your solution without
>>>>> running Power flow?
>>>>>
>>>>> I am using an algorithm that changes the amount of demand
>>>>> and generation at each bus (Pd,Pg) to mitigate over load of lines and
>>>>> enhance the security. The output of my algorithm (Pd,Pg) must be checked
>>>>> to
>>>>> insure that they are feasible. So, it is obvious that by assigning this
>>>>> new
>>>>> calculated values instead of old Pg and Pd (generation and demand ) at
>>>>> each
>>>>> bus and running a simple power flow (pf), I can find out if the calculated
>>>>> results are feasible.
>>>>> However, this is an optimization problem. Therefore, feasibility
>>>>> of the suggested values of Pd and Pg by fmincon must be verified at each
>>>>> iteration(1000 times). Basically, it can be accomplished by running a
>>>>> power
>>>>> flow in constraint function of fmincon and checking lines loading and
>>>>> voltage range.
>>>>> But, doing so is very time consuming and takes more than an
>>>>> hours to calculate the feasible Pg and Pd. So I was wonder if there is
>>>>> another faster way to check the feasibility of the calculated Pd,Pg
>>>>> without
>>>>> running power flow.
>>>>> I think, if the solution (calculated Pd,Pg) satisfies these
>>>>> constraints,the answer is then feasible and there is no need to run power
>>>>> flow.
>>>>> But I have some problems in considering the constraints 5,8,9,10.
>>>>>
>>>>> 1) Pgmin(j) <Pg(j) <Pgmax(j) % j indicates buses with generator
>>>>> 2) Qgmin(j) <Qg(j) <Qgmax(j) % j indicates buses with generator
>>>>> 3) Pdmin(k) <Pd(k) <Pdmax(k) % k indicates buses with demand
>>>>> 4) Qd(k)=Pd(k)*tan(phi(k)) % ph(k) is load P.F at bus K
>>>>> Nodal active balance:
>>>>> 5) Pg(n) -Pd(n)=|V(i)| sum( |Y(h,n)|
>>>>> .V(h)*Cos(del(n)-del(h)-teta(n,h)))
>>>>> where Pg(n) is nodal generation which is calculated as:
>>>>> 6) Pg(n)=sum (Pg(j)) at bus n
>>>>> 7) Pd(n)=sum (Pd(k)) at bus n
>>>>> teta(n,h) is the angle of Y(n,h), del(n) and del(h) are bust angle
>>>>> of n and h.
>>>>> Nodal reactive balance:
>>>>> 8) Qg(n) -Qd(n)=|V(i)| sum( |Y(h,n)|
>>>>> .V(h)*Sin(del(n)-del(h)-teta(n,h)))
>>>>>
>>>>> 9) Vmin(n)<V(n)<Vmax(n) % voltage range at all buses.
>>>>> 10) S(i,j)<Smax(i,j). transmitted apparent power from i,j must be
>>>>> less than capacity of line.
>>>>>
>>>>> Thanks in advance
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>
>>
>>
>
