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 <r...@cornell.edu> 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 <electricaltranslat...@gmail.com> > 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 <electricaltranslat...@gmail.com > > 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 <r...@cornell.edu> 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. <abhy...@mcs.anl.gov> >>> wrote: >>> >>> http://www.pserc.cornell.edu//matpower/docs/ref/ >>> >>> On Apr 9, 2015, at 4:00 PM, "Electric" <electricaltranslat...@gmail.com> >>> 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. < >>> abhy...@mcs.anl.gov> 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:* bounce-119032919-33970...@list.cornell.edu [ >>>> bounce-119032919-33970...@list.cornell.edu] on behalf of Electric [ >>>> electricaltranslat...@gmail.com] >>>> *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 >>>> >>>> >>>> >>> >>> >> > >