$ontext Filename : OPF-Formulations.gms Author : Steffen Rebennack & Stephen Frank Release : Sep 19, 2012 Version : 1.00 This program reads in an OPF instance. The electric network contains 5 buses, 6 branches and one phase-shifting and tap-changing transformer. The three different power flow models are established and solved with different solvers available in GAMS. This GAMS file is based on the following publication: S. Frank and S. Rebennack "An Introduction to Optimal Power Flow: Theory, Formulation, and Examples" IIE Transactions Notation and equation references are give with respect to this publication. --MODEL COMPARISON-------------------------------------------------------- | Function types macros taps | add./sub. mult. trig. inverse trig. OPF_PolarRecta no no | x x x yes no | x x x no yes | x x x yes yes | x x x OPF_PolarPolar no no | x x x yes no | x x x no yes | x x x x yes yes | x x x x OPF_RectaRecta no no | x x yes no | x x no yes | x x x yes yes | x x x --MODEL RESULTS----------------------------------------------------------- model OPF_PolarRecta results: with controllable taps without controllable taps macro no macro macro no macro SNOPT 0.4016596 0.4016595 0.4041438 0.4041438 MINOS 0.4016596 0.4016596 0.4041438 0.4041438 CONOPT 0.4016596 0.4016596 0.4041438 0.4041438 LINDOGlobal 0.4016596+ 0.4016596+ 0.4041438* 0.4041438* model OPF_PolarPolar results: with controllable taps without controllable taps macro no macro macro no macro SNOPT 0.4016596 0.4016596 0.4041438 0.4041438 MINOS 0.4016596 0.4016596 0.4041438 0.4041438 CONOPT 0.4016596 0.4016596 0.4041438 0.4041438 LINDOGlobal - - - 0.4041438* model OPF_RectaRecta results: with controllable taps without controllable taps macro no macro macro no macro SNOPT 0.4016596 0.4016596 0.4041438 0.4041438 MINOS 0.4016596 0.4016596 0.4041438 0.4041438 CONOPT 0.4016596 0.4016596 0.4041438 0.4041438 LINDOGlobal 0.4016596+ 0.4016596+ 0.4041438* 0.4041438* * proven global optimal solution + proven global optimal solution within optimality tolerance of 1% - solver reports that the problem is infeasible NOTES: 1. For the cases with controllable taps, LINDOGlobal finds the same optimum in the initial search as the other solvers, but requires an impractically long runtime to reduce the optimality gap to the default value of 1e-6. Therefore, for these model cases, an options file 'lindoglobal.opt' was created containing the following lines: * Set a reduced optimality tolerance OPTTOL 1e-2 This increases the default optimality tolerance to 1%, a level at which LINDOGlobal will terminate successfully. The 'lindoglobal.opt' file, if it exists, is automatically enabled when the compile-time switch %UseControllableTaps% is set to 'yes'. $offtext ** Compile-time controls * end-of-line comment symbol $eolcom // * Define whether we use macros for the admittance matrix elements: * yes (use macros), no (use auxiliary variables) $ set UseAdmittanceMacros yes * Define whether to activate transformer taps and phase angles * as decision variables: yes, no $ set UseControllableTaps yes * Define whether to turn on diagnostic displays of variables: yes, no $ set DiagDisplay yes * ============================================================================= * Data * ============================================================================= * ----------------------------------------------------------------------------- * Sets * ----------------------------------------------------------------------------- Set i set of system buses (nodes) /1*5/; * create copies Alias (i,k); * Display sets $IFTHEN %DiagDisplay%==yes display i; $ENDIF * ----------------------------------------------------------------------------- * Data * ----------------------------------------------------------------------------- Parameter g(i,k) 'conductance of branch (i,k)' b(i,k) 'susceptance of branch (i,k)'; * define the bus data: Table 5.1 Table BusData(i,*) Bus data follows PL QL VMin VMax gS bS 1 1.00 1.00 2 0.95 1.05 0.30 3 0.95 1.05 0.05 4 0.900 0.400 0.95 1.05 5 0.239 0.129 0.95 1.05 ; * define the generation data: Table 5.2 Table GenerationData(i,*) Generator data follows PGMin PGMax QGMin QGMax 1 -inf inf -inf inf 3 0.10 0.40 -0.20 0.30 4 0.05 0.40 -0.20 0.20 ; * define the branch data: Table 4.1 Table BranchData(i,k,*) Branch data follows R X b T phi 1.2 0.300 1.3 0.023 0.145 0.040 2.4 0.006 0.032 0.010 3.4 0.020 0.260 -3.0 3.5 0.320 0.98 4.5 0.500 ; // Define nominal tap ratio as 1.0 for all branches without a different // value specified. BranchData(i,k,'T')$(not BranchData(i,k,'T')) = 1; // Convert phase angle from degrees to radians BranchData(i,k,'phi')=(pi/180)*BranchData(i,k,'phi'); * define controllable phase-shifting transformer limits Table PhaseAngleLimits(i,k,*) Phase-shifting transformer phase angle limits phiMin phiMax 3.4 -30.0 30.0 ; // Convert phase angle limits from degrees to radians PhaseAngleLimits(i,k,'phiMin')=(pi/180)*PhaseAngleLimits(i,k,'phiMin'); PhaseAngleLimits(i,k,'phiMax')=(pi/180)*PhaseAngleLimits(i,k,'phiMax'); * define controllable tap-changing transformer limits Table TapRatioLimits(i,k,*) Tap-changing transformer voltage tap limits TMin TMax 3.5 0.95 1.05 ; * define the cost data Table C(i,*) linear quadratic 1 0.35 3 0.20 0.40 4 0.30 0.50 ; * Display data $IFTHEN %DiagDisplay%==yes display BusData, GenerationData, BranchData, PhaseAngleLimits, TapRatioLimits, C; $ENDIF * ----------------------------------------------------------------------------- * Derived Sets * ----------------------------------------------------------------------------- Set setG(i) set of controllable generation buses setL(i,k) set of system branches (arcs) setH(i,k) set of branches with controllable phase-shifting transformers setK(i,k) set of branches with controllable tap-changing transformers setNZ(i,k) set of non-zero admittance matrix entries setNF(i,k) set of variable (non-fixed) admittance matrix entries; * Compute set of lines setL(i,k)$(BranchData(i,k,'R') or BranchData(i,k,'X'))=1; * Compute set of generators setG(i)$(GenerationData(i,'PGMin'))=1; $IFTHEN %UseControllableTaps%==yes // Activate sets for controllable taps setH(i,k)$(PhaseAngleLimits(i,k,'phiMin') or PhaseAngleLimits(i,k,'phiMax'))=1; setK(i,k)$(TapRatioLimits(i,k,'TMin') or TapRatioLimits(i,k,'TMax'))=1; $ELSE // Use empty sets for controllable taps setH(i,k)=0; setK(i,k)=0; $ENDIF * Compute set of non-zero admittance matrix entries setNZ(i,i)$(sum(k,setL(i,k))+sum(k,setL(k,i)))=1; setNZ(i,k)$setL(i,k)=1; setNZ(i,k)$setL(k,i)=1; * Compute set of non-fixed admittance matrix entries setNF(i,i)$(sum(k,setH(i,k))+sum(k,setH(k,i)) +sum(k,setK(k,i))+sum(k,setK(k,i)))=1; setNF(i,k)$(setH(i,k) or setH(k,i))=1; setNF(i,k)$(setK(i,k) or setK(k,i))=1; * Display derived sets $IFTHEN %DiagDisplay%==yes display setL, setG, setH, setK, setNZ, setNF; $ENDIF * ----------------------------------------------------------------------------- * Derived Data * ----------------------------------------------------------------------------- * Calculate series admittance of each branch using equations (4.5) * Series conductance g(i,k)$setL(i,k)=BranchData(i,k,'R')/(BranchData(i,k,'R') *BranchData(i,k,'R')+BranchData(i,k,'X')*BranchData(i,k,'X')); * Series susceptance b(i,k)$setL(i,k)=-BranchData(i,k,'X')/(BranchData(i,k,'R') *BranchData(i,k,'R')+BranchData(i,k,'X')*BranchData(i,k,'X')); * Display derived data $IFTHEN %DiagDisplay%==yes display g, b; $ENDIF * ============================================================================= * Admittance Matrix Macros * ============================================================================= * Define macros to calculate admittance matrix elements $IFTHEN %UseAdmittanceMacros%==yes * Conductance matrix entries $ macro GbusD(i) ( BusData(i,'gS')+sum(k$setL(i,k),sqr(Tinv(i,k))*(g(i,k)))+sum(k$setL(k,i),g(k,i)) ) $ macro GbusO(i,k) ( -Tinv(i,k)*(g(i,k)*cos(phi(i,k))-b(i,k)*sin(phi(i,k)))-Tinv(k,i)*(g(k,i)*cos(phi(k,i))+b(k,i)*sin(phi(k,i))) ) * Susceptance matrix entries $ macro BbusD(i) ( BusData(i,'bS')+sum(k$setL(i,k),sqr(Tinv(i,k))*(b(i,k)+BranchData(i,k,'b')/2))+sum(k$setL(k,i),(b(k,i)+BranchData(k,i,'b')/2)) ) $ macro BbusO(i,k) ( -Tinv(i,k)*(g(i,k)*sin(phi(i,k))+b(i,k)*cos(phi(i,k)))-Tinv(k,i)*(-g(k,i)*sin(phi(k,i))+b(k,i)*cos(phi(k,i))) ) * Admittance matrix magnitude entries $ macro YbusD(i) ( sqrt( sqr(GbusD(i))+sqr(BbusD(i)) ) ) $ macro YbusO(i,k) ( sqrt( sqr(GbusO(i,k))+sqr(BbusO(i,k)) ) ) * Admittance matrix angle entries $ macro thetaD(i) ( arctan2(BbusD(i),GbusD(i)) ) $ macro thetaO(i,k) ( arctan2(BbusO(i,k),GbusO(i,k)) ) $ENDIF * ============================================================================= * Variables * ============================================================================= Variables * Voltages E(i) real component of complex voltage at bus i F(i) imaginary component of complex voltage at bus i delta(i) phase angle at bus i * Injected power P(i) net injected real power at bus i Q(i) net injected reactive power at bus i * Branch tap ratios phi(i,k) 'phase shift of complex tap ratio for branch (i,k)' Tinv(i,k) '1 / magnitude of complex tap ratio for branch (i,k)' * Aux. variables z objective $IFTHEN %UseAdmittanceMacros%==no Gbus(i,k) 'element (i,k) of the bus conductance matrix' Bbus(i,k) 'element (i,k) of the bus susceptance matrix' $ENDIF theta(i,k) 'angle of element (i,k) of the bus admittance matrix'; Positive variables * Voltages V(i) voltage magnitude * Aux. variables Ybus(i,k) 'magnitude of element (i,k) of the bus admittance matrix'; * ============================================================================== * Equations * ============================================================================== Equations z_obj objective function $IFTHEN %UseAdmittanceMacros%==no conductanceD(i) calculate diagonal entries of bus conductance matrix conductanceO(i,k) calculate off-diagonal entries of bus conductance matrix susceptanceD(i) calculate diagonal entries of bus susceptance matrix susceptanceO(i,k) calculate off-diagonal entries of bus susceptance matrix admitMag(i,k) constrain entries of bus admittance matrix (magnitude) admitAng(i,k) constrain entries of bus admittance matrix (angle) * admitAngSin(i,k) constrain entries of bus admittance matrix (angle) 1 of 2 * admitAngCos(i,k) constrain entries of bus admittance matrix (angle) 2 of 2 $ENDIF realPR(i) 'real power equations (polar voltage, rect. admittance)' reactivePR(i) 'reactive power equations (polar voltage, rect. admittance)' realPP(i) 'real power equations (polar voltage, polar admittance)' reactivePP(i) 'reactive power equations (polar voltage, polar admittance)' realRR(i) 'real power equations (rect. voltage, rect. admittance)' reactiveRR(i) 'reactive power equations (rect. voltage, rect. admittance)' voltminRR(i) 'lower voltage magnitude limit in rect. coordinates' voltmaxRR(i) 'upper voltage magnitude limit in rect. coordinates'; * objective function = cost of generated power z_obj.. z=e=sum(setG, C(setG,'linear')*(P(setG)+busData(setG,'PL')) + C(setG,'quadratic')*sqr(P(setG)+busData(setG,'PL')) ); $IFTHEN %UseAdmittanceMacros%==no * The following are required only if using auxiliary variables * (rather than macros) for admittance matrix element calculations. * Conductance matrix diagonal entries; equations (4.11) conductanceD(i)$setNF(i,i).. Gbus(i,i)=e=BusData(i,'gS') +sum(k$setL(i,k),sqr(Tinv(i,k))*(g(i,k))) +sum(k$setL(k,i),g(k,i)); * Conductance matrix off-diagonal entries; equations (4.12) conductanceO(i,k)$(setNF(i,k) and (ord(i)<>ord(k))).. Gbus(i,k)=e= -Tinv(i,k)*(g(i,k)*cos(phi(i,k))-b(i,k)*sin(phi(i,k))) -Tinv(k,i)*(g(k,i)*cos(phi(k,i))+b(k,i)*sin(phi(k,i))); * Susceptance matrix diagonal entries; equations (4.13) susceptanceD(i)$setNF(i,i).. Bbus(i,i)=e=BusData(i,'bS') +sum(k$setL(i,k),sqr(Tinv(i,k))*(b(i,k)+BranchData(i,k,'b')/2)) +sum(k$setL(k,i),(b(k,i)+BranchData(k,i,'b')/2)); * Susceptance matrix off-diagonal entries; equations (4.14) susceptanceO(i,k)$(setNF(i,k) and (ord(i)<>ord(k))).. Bbus(i,k)=e= -Tinv(i,k)*(g(i,k)*sin(phi(i,k))+b(i,k)*cos(phi(i,k))) -Tinv(k,i)*(-g(k,i)*sin(phi(k,i))+b(k,i)*cos(phi(k,i))); * Admittance matrix magnitudes admitMag(i,k)$setNF(i,k).. sqr(Ybus(i,k))=e=sqr(Gbus(i,k))+sqr(Bbus(i,k)); * Admittance matrix angles admitAng(i,k)$setNF(i,k).. theta(i,k)=e=arctan2(Bbus(i,k),Gbus(i,k)); *admitAngSin(i,k)$setNF(i,k).. * sin(theta(i,k))=e=Gbus(i,k)/Ybus(i,k); * *admitAngCos(i,k)$setNF(i,k).. * cos(theta(i,k))=e=Bbus(i,k)/Ybus(i,k); * The following are defined differently depending on whether * macros or aux. variables are used for the admittance matrix. ** Polar voltage coordinates, rectangular admittance coordinates * active power balance equations; equations (4.15) realPR(i).. P(i)=e=V(i)*sum(k,V(k)*(Gbus(i,k)*cos(delta(i)-delta(k)) +Bbus(i,k)*sin(delta(i)-delta(k)))); * reactive power balance equations; equations (4.16) reactivePR(i).. Q(i)=e=V(i)*sum(k,V(k)*(Gbus(i,k)*sin(delta(i)-delta(k)) -Bbus(i,k)*cos(delta(i)-delta(k)))); ** Polar voltage coordinates, polar admittance coordinates * active power balance equations; equations (4.17) realPP(i).. P(i)=e=V(i)*sum(k,V(k)*Ybus(i,k)*cos(delta(i)-delta(k)-theta(i,k))); * reactive power balance equations; equations (4.18) reactivePP(i).. Q(i)=e=V(i)*sum(k,V(k)*Ybus(i,k)*sin(delta(i)-delta(k)-theta(i,k))); ** Rectangular voltage coordinates, rectangular admittance coordinates * active power balance equations; equations (4.19) realRR(i).. P(i)=e=sum(k,Gbus(i,k)*(E(i)*E(k)+F(i)*F(k)) +Bbus(i,k)*(F(i)*E(k)-E(i)*F(k))); * reactive power balance equations; equations (4.20) reactiveRR(i).. Q(i)=e=sum(k,Gbus(i,k)*(F(i)*E(k)-E(i)*F(k)) -Bbus(i,k)*(E(i)*E(k)+F(i)*F(k))); $ELSE ** Polar voltage coordinates, rectangular admittance coordinates * active power balance equations; equations (4.15) realPR(i).. P(i)=e=V(i)*sum(k$(ord(i)<>ord(k)),V(k)*(GbusO(i,k)*cos(delta(i)-delta(k)) +BbusO(i,k)*sin(delta(i)-delta(k))))+sqr(V(i))*GbusD(i); * reactive power balance equations; equations (4.16) reactivePR(i).. Q(i)=e=V(i)*sum(k$(ord(i)<>ord(k)),V(k)*(GbusO(i,k)*sin(delta(i)-delta(k)) -BbusO(i,k)*cos(delta(i)-delta(k))))-sqr(V(i))*BbusD(i); ** Polar voltage coordinates, polar admittance coordinates * active power balance equations; equations (4.17) realPP(i).. P(i)=e=V(i)*sum(k$((ord(i)<>ord(k)) and setNZ(i,k)), V(k)*YbusO(i,k)*cos(delta(i)-delta(k)-thetaO(i,k))) +sqr(V(i))*YbusD(i)*cos(-thetaD(i)); * reactive power balance equations; equations (4.18) reactivePP(i).. Q(i)=e=V(i)*sum(k$((ord(i)<>ord(k)) and setNZ(i,k)), V(k)*YbusO(i,k)*sin(delta(i)-delta(k)-thetaO(i,k))) +sqr(V(i))*YbusD(i)*sin(-thetaD(i)); ** Rectangular voltage coordinates, rectangular admittance coordinates * active power balance equations; equations (4.19) realRR(i).. P(i)=e=sum(k$(ord(i)<>ord(k)),GbusO(i,k)*(E(i)*E(k)+F(i)*F(k)) +BbusO(i,k)*(F(i)*E(k)-E(i)*F(k))) +GbusD(i)*(E(i)*E(i)+F(i)*F(i)); * reactive power balance equations; equations (4.20) reactiveRR(i).. Q(i)=e=sum(k$(ord(i)<>ord(k)),GbusO(i,k)*(F(i)*E(k)-E(i)*F(k)) -BbusO(i,k)*(E(i)*E(k)+F(i)*F(k))) -BbusD(i)*(E(i)*E(i)+F(i)*F(i)); $ENDIF * Voltage magnitude limits (rectangular voltage coordinates) voltminRR(i).. E(i)*E(i)+F(i)*F(i)=g=sqr(BusData(i,'VMin')); voltmaxRR(i).. E(i)*E(i)+F(i)*F(i)=l=sqr(BusData(i,'VMax')); * ============================================================================= * Variable Bounds * ============================================================================= * Net real power (generation - load) P.lo(setG)=GenerationData(setG,'PGMin')-busData(setG,'PL'); P.up(setG)=GenerationData(setG,'PGMax')-busData(setG,'PL'); * Net reactive power (generation - load) Q.lo(setG)=GenerationData(setG,'QGMin')-busData(setG,'QL'); Q.up(setG)=GenerationData(setG,'QGMax')-busData(setG,'QL'); * Voltage magnitude (polar coordinates) V.lo(i)=BusData(i,'VMin'); V.up(i)=BusData(i,'VMax'); * Voltage magnitude (rectangular coordinates) E.lo(i)=-BusData(i,'VMax'); E.up(i)=BusData(i,'VMax'); F.lo(i)=-BusData(i,'VMax'); F.up(i)=BusData(i,'VMax'); * Voltage angle delta.lo(i)=-pi; delta.up(i)=+pi; * Branch phase angle phi.lo(i,k)$setH(i,k)=PhaseAngleLimits(i,k,'phiMin'); phi.up(i,k)$setH(i,k)=PhaseAngleLimits(i,k,'phiMax'); * 1 / branch tap ratio -- note: the lower limit of Tinv corresponds to 'TMax' Tinv.lo(i,k)$setK(i,k)=1/TapRatioLimits(i,k,'TMax'); Tinv.up(i,k)$setK(i,k)=1/TapRatioLimits(i,k,'TMin'); * Admittance matrix angle theta.lo(i,k)=-pi; theta.up(i,k)=+pi; * ============================================================================= * Fixed Variables * ============================================================================= * power at buses without generation P.fx(i)$(not setG(i))=-busData(i,'PL'); Q.fx(i)$(not setG(i))=-busData(i,'QL'); * non-controllable transformer voltage ratios Tinv.fx(i,k)$(not setK(i,k))=1/BranchData(i,k,'T'); * non-controllable transformer phase angles phi.fx(i,k)$(not setH(i,k))=BranchData(i,k,'phi'); * voltage angle at slack bus (bus 1) delta.fx('1')=0; F.fx('1')=0; $IFTHEN %UseAdmittanceMacros%==no * conductance matrix entries not influenced by controllable transformers Gbus.fx(i,i)$(not(setNF(i,i)))= BusData(i,'gS') +sum(k$setL(i,k),sqr(Tinv.l(i,k))*(g(i,k))) +sum(k$setL(k,i),g(k,i)); Gbus.fx(i,k)$(not(setNF(i,k)) and (ord(i)<>ord(k)))= -Tinv.l(i,k)*(g(i,k)*cos(phi.l(i,k))-b(i,k)*sin(phi.l(i,k))) -Tinv.l(k,i)*(g(k,i)*cos(phi.l(k,i))+b(k,i)*sin(phi.l(k,i))); * susceptance matrix entries not influenced by controllable transformers Bbus.fx(i,i)$(not(setNF(i,i)))= BusData(i,'bS') +sum(k$setL(i,k),sqr(Tinv.l(i,k))*(b(i,k)+BranchData(i,k,'b')/2)) +sum(k$setL(k,i),(b(k,i)+BranchData(k,i,'b')/2)); Bbus.fx(i,k)$(not(setNF(i,k)) and (ord(i)<>ord(k)))= -Tinv.l(i,k)*(g(i,k)*sin(phi.l(i,k))+b(i,k)*cos(phi.l(i,k))) -Tinv.l(k,i)*(-g(k,i)*sin(phi.l(k,i))+b(k,i)*cos(phi.l(k,i))); * admittance matrix entries in polar coordinates that are fixed to zero Ybus.fx(i,k)$(not(setNZ(i,k)))=0; theta.fx(i,k)$(not(setNZ(i,k)))=0; * admittance matrix magnitude entries not influenced by controllable transformers Ybus.fx(i,k)$(setNZ(i,k) and not(setNF(i,k)))= sqrt(sqr(Gbus.l(i,k))+sqr(Bbus.l(i,k))); * admittance matrix angle entries not influenced by controllable transformers theta.fx(i,k)$(setNZ(i,k) and not(setNF(i,k)))= arctan2(Bbus.l(i,k),Gbus.l(i,k)); $ENDIF * ============================================================================= * Initial Variable Values * ============================================================================= * voltage magnitude V.l(i)=(BusData(i,'VMin')+BusData(i,'VMax'))/2; E.l(i)=V.l(i); F.l(i)=0; * voltage angle delta.l(i)=0; * controllable transformer voltage ratios Tinv.l(i,k)$setK(i,k)=(Tinv.lo(i,k)+Tinv.up(i,k))/2; * controllable transformer phase angles phi.l(i,k)$setH(i,k)=(phi.lo(i,k)+phi.up(i,k))/2; * admittance matrix magnitudes (to avoid division by zero) Ybus.l(i,k)$setNF(i,k)=1; $IFTHEN %UseAdmittanceMacros%==no * conductance matrix entries influenced by controllable transformers Gbus.l(i,i)$setNF(i,i)= BusData(i,'gS') +sum(k$setL(i,k),sqr(Tinv.l(i,k))*(g(i,k))) +sum(k$setL(k,i),g(k,i)); Gbus.l(i,k)$(setNF(i,k) and (ord(i)<>ord(k)))= -Tinv.l(i,k)*(g(i,k)*cos(phi.l(i,k))-b(i,k)*sin(phi.l(i,k))) -Tinv.l(k,i)*(g(k,i)*cos(phi.l(k,i))+b(k,i)*sin(phi.l(k,i))); * susceptance matrix entries influenced by controllable transformers Bbus.l(i,i)$setNF(i,i)= BusData(i,'bS') +sum(k$setL(i,k),sqr(Tinv.l(i,k))*(b(i,k)+BranchData(i,k,'b')/2)) +sum(k$setL(k,i),(b(k,i)+BranchData(k,i,'b')/2)); Bbus.l(i,k)$(setNF(i,k) and (ord(i)<>ord(k)))= -Tinv.l(i,k)*(g(i,k)*sin(phi.l(i,k))+b(i,k)*cos(phi.l(i,k))) -Tinv.l(k,i)*(-g(k,i)*sin(phi.l(k,i))+b(k,i)*cos(phi.l(k,i))); * admittance matrix magnitude entries influenced by controllable transformers Ybus.l(i,k)$setNF(i,k)= sqrt(sqr(Gbus.l(i,k))+sqr(Bbus.l(i,k))); * admittance matrix angle entries influenced by controllable transformers theta.l(i,k)$setNF(i,k)= arctan2(Bbus.l(i,k),Gbus.l(i,k)); $ENDIF * Display initial variable values $IFTHEN %DiagDisplay%==yes $ontext Starting values of fixed and free variables follow: $offtext display V.l, delta.l, E.l, F.l, P.l, Q.l, Tinv.l, phi.l; $IFTHEN.TWO %UseAdmittanceMacros%==no display Gbus.l, Bbus.l, Ybus.l, theta.l; $ENDIF.TWO $ENDIF * ============================================================================= * Model Definitions * ============================================================================= * The constraints included in the model depend on whether or not aux. variables * and their corresponding equations are used for the admittance matrix. $IFTHEN %UseAdmittanceMacros%==no * polar coordinates for voltage, rectangular coordinates for admittance Model OPF_PolarRecta /z_obj,conductanceD,conductanceO,susceptanceD, susceptanceO,realPR,reactivePR/; * polar coordinates for voltage, polar coordinates for admittance Model OPF_PolarPolar /z_obj,conductanceD,conductanceO,susceptanceD, susceptanceO,admitMag,admitAng, realPP,reactivePP/; * rectangular coordinates for voltage, rectangular coordinates for admittance Model OPF_RectaRecta /z_obj,conductanceD,conductanceO,susceptanceD, susceptanceO,realRR,reactiveRR, voltminRR,voltmaxRR/; $ELSE * polar coordinates for voltage, rectangular coordinates for admittance Model OPF_PolarRecta /z_obj,realPR,reactivePR/; * polar coordinates for voltage, polar coordinates for admittance Model OPF_PolarPolar /z_obj,realPP,reactivePP/; * rectangular coordinates for voltage, rectangular coordinates for admittance Model OPF_RectaRecta /z_obj,realRR,reactiveRR,voltminRR,voltmaxRR/; $ENDIF * The models, ranked in order of increasing nonlinearity: * 1. OPF_RectaRecta * functions: addition, multiplication, trigonometric (only if taps included) * 2. OPF_PolarRecta * functions: addition, multiplication, trigonometric * 3. OPF_PolarPolar * terms: linear, multiplication, division, trigonometric, * inverse trigonometric (only if taps included) * ============================================================================= * Solver Settings * ============================================================================= * set solver options option limrow = 1e9 ; // show all constraints in the list file option limcol = 0 ; // no display of columns in the *.lst file option optcr = 0.00001 ; // global solver needs small relative gap or zero option reslim = 9999999 ; * Loads solver options from external file; used to reduce required optimality * tolerance for LINDOGlobal for cases with controllable taps. $IFTHEN %UseControllableTaps%==yes OPF_PolarRecta.optfile = 1 ; OPF_PolarPolar.optfile = 1 ; OPF_RectaRecta.optfile = 1 ; $ENDIF * select a global solver (as we have a non-convex objective function) * option NLP=SNOPT ; // solves to global optimum * option NLP=MINOS ; // solves to global optimum * option NLP=CONOPT ; // solves to global optimum option NLP=LINDOGlobal ; // solves to global optimum * ============================================================================= * Solve Commands * ============================================================================= * solve the model * Solve OPF_PolarRecta using NLP min z ; * Solve OPF_PolarPolar using NLP min z ; Solve OPF_RectaRecta using NLP min z ; ** Calculate aux. outputs * Calculate the voltage angle output in degrees Parameter deltaDegree(i) ; deltaDegree(i) = delta.l(i) / pi * 180 ; * Calculate the tap ratios as originally specified Parameter T(i,k) ; T(i,k)$setL(i,k) = 1 / Tinv.l(i,k) ; * Calculate the branch phase angles in degrees Parameter phiDegree(i,k) ; phiDegree(i,k)$setL(i,k) = phi.l(i,k) / pi * 180 ; * Calculate the generator powers Parameters PG(i), QG(i); PG(i) = P.l(i) + busData(i,'PL'); QG(i) = Q.l(i) + busData(i,'QL') * Display aux. outputs $IFTHEN %DiagDisplay%==yes display deltaDegree, T, phiDegree, PG, QG; $ENDIF * write the solution into a file *execute_unload "OPF-results.gdx" Bbus.l, Gbus.l, P.l, Q.l, execute_unload "OPF-results.gdx" P.l, Q.l, V.l, delta.l, T, phiDegree, deltaDegree ; * =================== end of OPF-Formulations.gms =============================