| 1 | package BugWithHomotopy
|
|---|
| 2 |
|
|---|
| 3 | package Examples
|
|---|
| 4 |
|
|---|
| 5 | model NetworkModelWithHomotopy
|
|---|
| 6 | extends Modelica.Icons.Example;
|
|---|
| 7 | extends BugWithHomotopy.Examples.NetworkModelWithoutHomotopy(redeclare Elements.CompositeLoadWithDummyHomotopy composite_load);
|
|---|
| 8 | annotation(
|
|---|
| 9 | experiment(StartTime = 0, StopTime = 10, Tolerance = 1e-06, Interval = 0.01),
|
|---|
| 10 | __OpenModelica_commandLineOptions = "--matchingAlgorithm=PFPlusExt --indexReductionMethod=dynamicStateSelection -d=initialization,NLSanalyticJacobian,newInst",
|
|---|
| 11 | __OpenModelica_simulationFlags(lv = "LOG_NLS,LOG_NLS_V,LOG_STATS", s = "ida"));
|
|---|
| 12 | end NetworkModelWithHomotopy;
|
|---|
| 13 | model NetworkModelWithoutHomotopy
|
|---|
| 14 | extends Modelica.Icons.Example;
|
|---|
| 15 | inner BugWithHomotopy.Elements.SystemBase SysData(SBase = 1e+6, fBase = 50) annotation(
|
|---|
| 16 | Placement(visible = true, transformation(origin = {-70, 90}, extent = {{-30, -10}, {30, 10}}, rotation = 0)));
|
|---|
| 17 |
|
|---|
| 18 | parameter Real P_0 = 18.549341834292655e3;
|
|---|
| 19 | parameter Real Q_0 = 6.258853e3;
|
|---|
| 20 | parameter Real v_0 = 1.0003343 * 20e3;
|
|---|
| 21 | parameter Real angle_0 = 1;
|
|---|
| 22 | BugWithHomotopy.Elements.ElmVac source(VBase = 20e3, angle_0 = angle_0, v_0 = v_0, P_0 = -P_0, Q_0 = -Q_0) annotation(
|
|---|
| 23 | Placement(visible = true, transformation(origin = {-36, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
|
|---|
| 24 | Modelica.Blocks.Sources.Constant synchronous_frequency(k = 1) annotation(
|
|---|
| 25 | Placement(visible = true, transformation(origin = {-70, 30}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
|
|---|
| 26 | BugWithHomotopy.Elements.CompositeLoad composite_load(P_0=P_0, Q_0=Q_0, VBase = 20e3, angle_0=angle_0, omega_0 = 1, v_0=v_0) annotation(
|
|---|
| 27 | Placement(visible = true, transformation(origin = {0, 0}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
|
|---|
| 28 | Modelica.Blocks.Sources.Step voltage(height = -0.01, offset = v_0 / 20e3, startTime = 1) annotation(
|
|---|
| 29 | Placement(visible = true, transformation(origin = {-70, -30}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
|
|---|
| 30 | equation
|
|---|
| 31 | connect(source.p, composite_load.p) annotation(
|
|---|
| 32 | Line(points = {{-36, 0}, {0, 0}}, color = {0, 85, 0}));
|
|---|
| 33 | connect(composite_load.win, synchronous_frequency.y) annotation(
|
|---|
| 34 | Line(points = {{10, 0}, {20, 0}, {20, 30}, {-58, 30}}, color = {0, 0, 127}));
|
|---|
| 35 | connect(source.fPu, synchronous_frequency.y) annotation(
|
|---|
| 36 | Line(points = {{-46, 6}, {-52, 6}, {-52, 30}, {-58, 30}}, color = {0, 0, 127}));
|
|---|
| 37 | connect(voltage.y, source.vPu) annotation(
|
|---|
| 38 | Line(points = {{-58, -30}, {-52, -30}, {-52, -4}, {-46, -4}}, color = {0, 0, 127}));
|
|---|
| 39 | annotation(
|
|---|
| 40 | experiment(StartTime = 0, StopTime = 10, Tolerance = 1e-06, Interval = 0.01),
|
|---|
| 41 | __OpenModelica_commandLineOptions = "--matchingAlgorithm=PFPlusExt --indexReductionMethod=dynamicStateSelection -d=initialization,NLSanalyticJacobian,newInst",
|
|---|
| 42 | __OpenModelica_simulationFlags(lv = "LOG_NLS,LOG_NLS_V,LOG_STATS", s = "ida"));
|
|---|
| 43 | end NetworkModelWithoutHomotopy;
|
|---|
| 44 | extends Modelica.Icons.ExamplesPackage;
|
|---|
| 45 |
|
|---|
| 46 | end Examples;
|
|---|
| 47 | extends Modelica.Icons.Package;
|
|---|
| 48 | import SI = Modelica.SIunits;
|
|---|
| 49 | import C = Modelica.Constants;
|
|---|
| 50 | import CM = Modelica.ComplexMath;
|
|---|
| 51 |
|
|---|
| 52 | package Elements
|
|---|
| 53 |
|
|---|
| 54 | model ElmVac
|
|---|
| 55 | extends Elements.BaseClasses.OnePort;
|
|---|
| 56 | SI.Angle phiu;
|
|---|
| 57 | Modelica.Blocks.Interfaces.RealInput fPu annotation(
|
|---|
| 58 | Placement(visible = true, transformation(origin = {-100, 50}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-100, 50}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
|
|---|
| 59 | Modelica.Blocks.Interfaces.RealInput vPu annotation(
|
|---|
| 60 | Placement(visible = true, transformation(origin = {-100, -50}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-100, -50}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
|
|---|
| 61 | initial equation
|
|---|
| 62 | phiu = angle_0;
|
|---|
| 63 | equation
|
|---|
| 64 | der(phiu) = 2 * fBase * C.pi * (fPu - 1);
|
|---|
| 65 | p.v = VBase * Complex(vPu * cos(phiu), vPu * sin(phiu));
|
|---|
| 66 | annotation(
|
|---|
| 67 | Icon(graphics = {Ellipse(fillColor = {255, 255, 255}, fillPattern = FillPattern.Solid, extent = {{-100, 100}, {100, -100}}, endAngle = 360), Line(origin = {-5.53254, -0.428994}, points = {{-40, 0}, {-20, 20}, {20, -20}, {40, 0}, {40, 0}})}));
|
|---|
| 68 | end ElmVac;
|
|---|
| 69 | model CompositeLoad
|
|---|
| 70 | extends BugWithHomotopy.Elements.BaseClasses.OnePort(p(i(re(nominal = M_b / VBase), im(nominal = M_b / VBase))));
|
|---|
| 71 | parameter SI.PerUnit omega_0 = 1 "Initial frequency" annotation(
|
|---|
| 72 | Dialog(group = "Power flow data"));
|
|---|
| 73 | // General Parameters
|
|---|
| 74 | parameter Real Kpm = 0.8 annotation(
|
|---|
| 75 | Dialog(group = "General Parameters"));
|
|---|
| 76 | parameter Real Mlf = 0.8 annotation(
|
|---|
| 77 | Dialog(group = "General Parameters"));
|
|---|
| 78 | // Parmeters of the exp model
|
|---|
| 79 | parameter Real expP = 2 annotation(
|
|---|
| 80 | Dialog(group = "Exponential Load Model Parameters"));
|
|---|
| 81 | parameter Real expQ = 2 annotation(
|
|---|
| 82 | Dialog(group = "Exponential Load Model Parameters"));
|
|---|
| 83 | // Parameters of the induction machine
|
|---|
| 84 | parameter SI.PerUnit Rs = 0.028 annotation(
|
|---|
| 85 | Dialog(group = "Machine Parameters"));
|
|---|
| 86 | parameter SI.PerUnit Xs = 0.01 annotation(
|
|---|
| 87 | Dialog(group = "Machine Parameters"));
|
|---|
| 88 | parameter SI.PerUnit Xm = 4.236 annotation(
|
|---|
| 89 | Dialog(group = "Machine Parameters"));
|
|---|
| 90 | parameter SI.PerUnit Rr = 0.00489 annotation(
|
|---|
| 91 | Dialog(group = "Machine Parameters"));
|
|---|
| 92 | parameter SI.PerUnit Xr = 0.1822 annotation(
|
|---|
| 93 | Dialog(group = "Machine Parameters"));
|
|---|
| 94 | parameter Real H = 1.67 annotation(
|
|---|
| 95 | Dialog(group = "Machine Parameters"));
|
|---|
| 96 | parameter Real expT = 2 annotation(
|
|---|
| 97 | Dialog(group = "Machine Parameters"));
|
|---|
| 98 | final parameter Real omega_base = 2 * C.pi * fBase;
|
|---|
| 99 | // Active and reactive powers
|
|---|
| 100 | final parameter Types.ActivePower Pm_0 = Kpm * P_0;
|
|---|
| 101 | final parameter Types.ApparentPower M_b = Pm_0 / Mlf;
|
|---|
| 102 | final parameter Types.ActivePower Pexp_0 = P_0 - Pm_0;
|
|---|
| 103 | final parameter Types.ReactivePower Qexp_0(fixed = false);
|
|---|
| 104 | final parameter Real x_ss = Xs + Xm;
|
|---|
| 105 | final parameter Real x_sr = Xm;
|
|---|
| 106 | final parameter Real x_rr = Xr + Xm;
|
|---|
| 107 | final parameter Real r_r = Rr;
|
|---|
| 108 | final parameter Real r_s = Rs;
|
|---|
| 109 | final parameter Real xhat = x_ss - x_sr / x_rr * x_sr;
|
|---|
| 110 | // Start parameters
|
|---|
| 111 | parameter Types.ComplexCurrent I_m0 = Complex(i_sd0, i_sq0) * M_b / VBase;
|
|---|
| 112 | parameter Types.ComplexCurrent I_e0 = I_0 - I_m0;
|
|---|
| 113 | final parameter SI.PerUnit u_sd0 = v_0 * cos(angle_0) / VBase;
|
|---|
| 114 | final parameter SI.PerUnit u_sq0 = v_0 * sin(angle_0) / VBase;
|
|---|
| 115 | final parameter SI.PerUnit psi_rd0(fixed = false);
|
|---|
| 116 | final parameter SI.PerUnit psi_rq0(fixed = false);
|
|---|
| 117 | final parameter SI.PerUnit psi_sq0(fixed = false);
|
|---|
| 118 | final parameter SI.PerUnit psi_sd0(fixed = false);
|
|---|
| 119 | final parameter SI.PerUnit i_sd0(start = I_0.re / (M_b / VBase), fixed = false);
|
|---|
| 120 | final parameter SI.PerUnit i_sq0(start = I_0.im / (M_b / VBase), fixed = false);
|
|---|
| 121 | final parameter SI.PerUnit w_0(start = omega_0, fixed = false);
|
|---|
| 122 | // Model variables
|
|---|
| 123 | // Angle w.r.t. the reference angle in the 50Hz rotating frame
|
|---|
| 124 | SI.PerUnit delta;
|
|---|
| 125 | // Synchronous frequency
|
|---|
| 126 | Modelica.Blocks.Interfaces.RealInput win annotation(
|
|---|
| 127 | Placement(visible = true, transformation(origin = {-100, 70}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-100, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
|
|---|
| 128 | Types.ComplexCurrent imach(re(start = I_m0.re), im(start = I_m0.im));
|
|---|
| 129 | Types.ComplexCurrent iexp(re(start = I_e0.re), im(start = I_e0.im));
|
|---|
| 130 | SI.PerUnit wr(start = w_0);
|
|---|
| 131 | SI.PerUnit psi_rd(start = psi_rd0);
|
|---|
| 132 | SI.PerUnit psi_rq(start = psi_rq0);
|
|---|
| 133 | SI.PerUnit psi_sd(start = psi_sd0);
|
|---|
| 134 | SI.PerUnit psi_sq(start = psi_sq0);
|
|---|
| 135 | SI.PerUnit i_sd(start = i_sd0);
|
|---|
| 136 | SI.PerUnit i_sq(start = i_sq0);
|
|---|
| 137 | SI.PerUnit u_sd(start = u_sd0);
|
|---|
| 138 | SI.PerUnit u_sq(start = u_sq0);
|
|---|
| 139 | SI.PerUnit T_m(start = (-i_sd0 * psi_sq0) + i_sq0 * psi_sd0);
|
|---|
| 140 | SI.PerUnit T_e(start = (-i_sd0 * psi_sq0) + i_sq0 * psi_sd0);
|
|---|
| 141 | SI.PerUnit T_0(start = (-i_sd0 * psi_sq0) + i_sq0 * psi_sd0);
|
|---|
| 142 | initial equation
|
|---|
| 143 | // Start values for the machine model
|
|---|
| 144 | 0 = r_r / x_rr * x_sr * i_sd0 - r_r / x_rr * psi_rd0 + (omega_0 - w_0) * psi_rq0;
|
|---|
| 145 | 0 = r_r / x_rr * x_sr * i_sq0 - r_r / x_rr * psi_rq0 - (omega_0 - w_0) * psi_rd0;
|
|---|
| 146 | psi_sd0 = xhat * i_sd0 + x_sr / x_rr * psi_rd0;
|
|---|
| 147 | psi_sq0 = xhat * i_sq0 + x_sr / x_rr * psi_rq0;
|
|---|
| 148 | u_sd0 = r_s * i_sd0 - omega_0 / 1 * psi_sq0;
|
|---|
| 149 | u_sq0 = r_s * i_sq0 + omega_0 / 1 * psi_sd0;
|
|---|
| 150 | Pm_0 / M_b = u_sd0 * i_sd0 + u_sq0 * i_sq0;
|
|---|
| 151 | // Initialization
|
|---|
| 152 | Qexp_0 = Q_0 - (u_sq * i_sd - u_sd * i_sq) * M_b;
|
|---|
| 153 | P = P_0;
|
|---|
| 154 | delta = 0;
|
|---|
| 155 | der(wr) = 0;
|
|---|
| 156 | der(psi_rd) = 0;
|
|---|
| 157 | der(psi_rq) = 0;
|
|---|
| 158 | equation
|
|---|
| 159 | // Integration of frequency differences to determine delta for d-q transformation
|
|---|
| 160 | der(delta) = omega_base * (win - 1);
|
|---|
| 161 | // d-q transformations
|
|---|
| 162 | [u_sd; u_sq] * VBase = [cos(delta), sin(delta); -sin(delta), cos(delta)] * [p.v.re; p.v.im];
|
|---|
| 163 | [i_sd; i_sq] * (M_b / VBase) = [cos(delta), sin(delta); -sin(delta), cos(delta)] * [imach.re; imach.im];
|
|---|
| 164 | // Mechanical equations
|
|---|
| 165 | der(T_0) = 0;
|
|---|
| 166 | T_m = T_0 * wr ^ expT;
|
|---|
| 167 | T_e = (-i_sd * psi_sq) + i_sq * psi_sd;
|
|---|
| 168 | der(wr) = 1 / (2 * H) * (T_e - T_m);
|
|---|
| 169 | // Electrical equations
|
|---|
| 170 | u_sd = r_s * i_sd - win / 1 * psi_sq;
|
|---|
| 171 | u_sq = r_s * i_sq + win / 1 * psi_sd;
|
|---|
| 172 | psi_sd = xhat * i_sd + x_sr / x_rr * psi_rd;
|
|---|
| 173 | psi_sq = xhat * i_sq + x_sr / x_rr * psi_rq;
|
|---|
| 174 | der(psi_rd) / omega_base = r_r / x_rr * x_sr * i_sd - r_r / x_rr * psi_rd + (win - wr) * psi_rq;
|
|---|
| 175 | der(psi_rq) / omega_base = r_r / x_rr * x_sr * i_sq - r_r / x_rr * psi_rq - (win - wr) * psi_rd;
|
|---|
| 176 | // Exponential load model
|
|---|
| 177 | p.v * CM.conj(iexp) = Complex(Pexp_0 * (CM.'abs'(p.v) / v_0) ^ expP, Qexp_0 * (CM.'abs'(p.v) / v_0) ^ expQ);
|
|---|
| 178 | p.i = imach + iexp;
|
|---|
| 179 | annotation(
|
|---|
| 180 | Icon(graphics = {Ellipse(fillColor = {255, 255, 255}, fillPattern = FillPattern.Solid, extent = {{-100, 100}, {100, -100}}, endAngle = 360), Text(origin = {-40.989, -17}, extent = {{-39.011, -3}, {118.989, -63}}, textString = "COMPOSITE")}, coordinateSystem(initialScale = 0.1)));
|
|---|
| 181 | end CompositeLoad;
|
|---|
| 182 |
|
|---|
| 183 | model CompositeLoadWithDummyHomotopy
|
|---|
| 184 | extends BugWithHomotopy.Elements.BaseClasses.OnePort(p(i(re(nominal = M_b / VBase), im(nominal = M_b / VBase))));
|
|---|
| 185 | parameter SI.PerUnit omega_0 = 1 "Initial frequency" annotation(
|
|---|
| 186 | Dialog(group = "Power flow data"));
|
|---|
| 187 | // General Parameters
|
|---|
| 188 | parameter Real Kpm = 0.8 annotation(
|
|---|
| 189 | Dialog(group = "General Parameters"));
|
|---|
| 190 | parameter Real Mlf = 0.8 annotation(
|
|---|
| 191 | Dialog(group = "General Parameters"));
|
|---|
| 192 | // Parmeters of the exp model
|
|---|
| 193 | parameter Real expP = 2 annotation(
|
|---|
| 194 | Dialog(group = "Exponential Load Model Parameters"));
|
|---|
| 195 | parameter Real expQ = 2 annotation(
|
|---|
| 196 | Dialog(group = "Exponential Load Model Parameters"));
|
|---|
| 197 | // Parameters of the induction machine
|
|---|
| 198 | parameter SI.PerUnit Rs = 0.028 annotation(
|
|---|
| 199 | Dialog(group = "Machine Parameters"));
|
|---|
| 200 | parameter SI.PerUnit Xs = 0.01 annotation(
|
|---|
| 201 | Dialog(group = "Machine Parameters"));
|
|---|
| 202 | parameter SI.PerUnit Xm = 4.236 annotation(
|
|---|
| 203 | Dialog(group = "Machine Parameters"));
|
|---|
| 204 | parameter SI.PerUnit Rr = 0.00489 annotation(
|
|---|
| 205 | Dialog(group = "Machine Parameters"));
|
|---|
| 206 | parameter SI.PerUnit Xr = 0.1822 annotation(
|
|---|
| 207 | Dialog(group = "Machine Parameters"));
|
|---|
| 208 | parameter Real H = 1.67 annotation(
|
|---|
| 209 | Dialog(group = "Machine Parameters"));
|
|---|
| 210 | parameter Real expT = 2 annotation(
|
|---|
| 211 | Dialog(group = "Machine Parameters"));
|
|---|
| 212 | final parameter Real omega_base = 2 * C.pi * fBase;
|
|---|
| 213 | // Active and reactive powers
|
|---|
| 214 | final parameter Types.ActivePower Pm_0 = Kpm * P_0;
|
|---|
| 215 | final parameter Types.ApparentPower M_b = Pm_0 / Mlf;
|
|---|
| 216 | final parameter Types.ActivePower Pexp_0 = P_0 - Pm_0;
|
|---|
| 217 | final parameter Types.ReactivePower Qexp_0(fixed = false);
|
|---|
| 218 | final parameter Real x_ss = Xs + Xm;
|
|---|
| 219 | final parameter Real x_sr = Xm;
|
|---|
| 220 | final parameter Real x_rr = Xr + Xm;
|
|---|
| 221 | final parameter Real r_r = Rr;
|
|---|
| 222 | final parameter Real r_s = Rs;
|
|---|
| 223 | final parameter Real xhat = x_ss - x_sr / x_rr * x_sr;
|
|---|
| 224 | // Start parameters
|
|---|
| 225 | parameter Types.ComplexCurrent I_m0 = Complex(i_sd0, i_sq0) * M_b / VBase;
|
|---|
| 226 | parameter Types.ComplexCurrent I_e0 = I_0 - I_m0;
|
|---|
| 227 | final parameter SI.PerUnit u_sd0 = v_0 * cos(angle_0) / VBase;
|
|---|
| 228 | final parameter SI.PerUnit u_sq0 = v_0 * sin(angle_0) / VBase;
|
|---|
| 229 | final parameter SI.PerUnit psi_rd0(fixed = false);
|
|---|
| 230 | final parameter SI.PerUnit psi_rq0(fixed = false);
|
|---|
| 231 | final parameter SI.PerUnit psi_sq0(fixed = false);
|
|---|
| 232 | final parameter SI.PerUnit psi_sd0(fixed = false);
|
|---|
| 233 | final parameter SI.PerUnit i_sd0(start = I_0.re / (M_b / VBase), fixed = false);
|
|---|
| 234 | final parameter SI.PerUnit i_sq0(start = I_0.im / (M_b / VBase), fixed = false);
|
|---|
| 235 | final parameter SI.PerUnit w_0(start = omega_0, fixed = false);
|
|---|
| 236 | // Model variables
|
|---|
| 237 | // Angle w.r.t. the reference angle in the 50Hz rotating frame
|
|---|
| 238 | SI.PerUnit delta;
|
|---|
| 239 | // Synchronous frequency
|
|---|
| 240 | Modelica.Blocks.Interfaces.RealInput win annotation(
|
|---|
| 241 | Placement(visible = true, transformation(origin = {-100, 70}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-100, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
|
|---|
| 242 | Types.ComplexCurrent imach(re(start = I_m0.re), im(start = I_m0.im));
|
|---|
| 243 | Types.ComplexCurrent iexp(re(start = I_e0.re), im(start = I_e0.im));
|
|---|
| 244 | SI.PerUnit wr(start = w_0);
|
|---|
| 245 | SI.PerUnit psi_rd(start = psi_rd0);
|
|---|
| 246 | SI.PerUnit psi_rq(start = psi_rq0);
|
|---|
| 247 | SI.PerUnit psi_sd(start = psi_sd0);
|
|---|
| 248 | SI.PerUnit psi_sq(start = psi_sq0);
|
|---|
| 249 | SI.PerUnit i_sd(start = i_sd0);
|
|---|
| 250 | SI.PerUnit i_sq(start = i_sq0);
|
|---|
| 251 | SI.PerUnit u_sd(start = u_sd0);
|
|---|
| 252 | SI.PerUnit u_sq(start = u_sq0);
|
|---|
| 253 | SI.PerUnit T_m(start = (-i_sd0 * psi_sq0) + i_sq0 * psi_sd0);
|
|---|
| 254 | SI.PerUnit T_e(start = (-i_sd0 * psi_sq0) + i_sq0 * psi_sd0);
|
|---|
| 255 | SI.PerUnit T_0(start = (-i_sd0 * psi_sq0) + i_sq0 * psi_sd0);
|
|---|
| 256 | initial equation
|
|---|
| 257 | // Start values for the machine model
|
|---|
| 258 | 0 = r_r / x_rr * x_sr * i_sd0 - r_r / x_rr * psi_rd0 + (omega_0 - w_0) * psi_rq0;
|
|---|
| 259 | 0 = r_r / x_rr * x_sr * i_sq0 - r_r / x_rr * psi_rq0 - (omega_0 - w_0) * psi_rd0;
|
|---|
| 260 | psi_sd0 = xhat * i_sd0 + x_sr / x_rr * psi_rd0;
|
|---|
| 261 | psi_sq0 = xhat * i_sq0 + x_sr / x_rr * psi_rq0;
|
|---|
| 262 | u_sd0 = r_s * i_sd0 - omega_0 / 1 * psi_sq0;
|
|---|
| 263 | u_sq0 = r_s * i_sq0 + omega_0 / 1 * psi_sd0;
|
|---|
| 264 | Pm_0 / M_b = u_sd0 * i_sd0 + u_sq0 * i_sq0;
|
|---|
| 265 | // Initialization
|
|---|
| 266 | Qexp_0 = Q_0 - (u_sq * i_sd - u_sd * i_sq) * M_b;
|
|---|
| 267 | P = P_0;
|
|---|
| 268 | delta = 0;
|
|---|
| 269 | der(wr) = 0;
|
|---|
| 270 | der(psi_rd) = 0;
|
|---|
| 271 | der(psi_rq) = 0;
|
|---|
| 272 | equation
|
|---|
| 273 | // Integration of frequency differences to determine delta for d-q transformation
|
|---|
| 274 | der(delta) = omega_base * (win - 1);
|
|---|
| 275 | // d-q transformations
|
|---|
| 276 | [u_sd; u_sq] * VBase = [cos(delta), sin(delta); -sin(delta), cos(delta)] * [p.v.re; p.v.im];
|
|---|
| 277 | [i_sd; i_sq] * (M_b / VBase) = [cos(delta), sin(delta); -sin(delta), cos(delta)] * [imach.re; imach.im];
|
|---|
| 278 | // Mechanical equations
|
|---|
| 279 | der(T_0) = 0;
|
|---|
| 280 | T_m = homotopy(T_0 * wr ^ expT, T_0 * wr ^ 0);
|
|---|
| 281 | T_e = (-i_sd * psi_sq) + i_sq * psi_sd;
|
|---|
| 282 | der(wr) = 1 / (2 * H) * (T_e - T_m);
|
|---|
| 283 | // Electrical equations
|
|---|
| 284 | u_sd = r_s * i_sd - win / 1 * psi_sq;
|
|---|
| 285 | u_sq = r_s * i_sq + win / 1 * psi_sd;
|
|---|
| 286 | psi_sd = xhat * i_sd + x_sr / x_rr * psi_rd;
|
|---|
| 287 | psi_sq = xhat * i_sq + x_sr / x_rr * psi_rq;
|
|---|
| 288 | der(psi_rd) / omega_base = r_r / x_rr * x_sr * i_sd - r_r / x_rr * psi_rd + (win - wr) * psi_rq;
|
|---|
| 289 | der(psi_rq) / omega_base = r_r / x_rr * x_sr * i_sq - r_r / x_rr * psi_rq - (win - wr) * psi_rd;
|
|---|
| 290 | // Exponential load model
|
|---|
| 291 | p.v * CM.conj(iexp) = Complex(Pexp_0 * (CM.'abs'(p.v) / v_0) ^ expP, Qexp_0 * (CM.'abs'(p.v) / v_0) ^ expQ);
|
|---|
| 292 | p.i = imach + iexp;
|
|---|
| 293 | annotation(
|
|---|
| 294 | Icon(graphics = {Ellipse(fillColor = {255, 255, 255}, fillPattern = FillPattern.Solid, extent = {{-100, 100}, {100, -100}}, endAngle = 360), Text(origin = {-40.989, -17}, extent = {{-39.011, -3}, {118.989, -63}}, textString = "COMPOSITE")}, coordinateSystem(initialScale = 0.1)));
|
|---|
| 295 | end CompositeLoadWithDummyHomotopy;
|
|---|
| 296 |
|
|---|
| 297 | record SystemBase "System Base Definition"
|
|---|
| 298 | parameter SI.ApparentPower SBase(displayUnit = "MVA") = 100e6 "System base" annotation(
|
|---|
| 299 | Evaluate = true);
|
|---|
| 300 | parameter SI.Frequency fBase = 50 "System frequency" annotation(
|
|---|
| 301 | Evaluate = true);
|
|---|
| 302 | annotation(
|
|---|
| 303 | Icon(coordinateSystem(preserveAspectRatio = false, extent = {{-120, -100}, {120, 100}}), graphics = {Rectangle(extent = {{-120, 100}, {120, -100}}, lineColor = {28, 108, 200}), Text(extent = {{-100, 40}, {100, 0}}, lineColor = {28, 108, 200}, horizontalAlignment = TextAlignment.Left, textString = "System Base: %S_b"), Text(extent = {{-100, -20}, {100, -60}}, lineColor = {28, 108, 200}, horizontalAlignment = TextAlignment.Left, textString = "Frequency: %fn"), Text(extent = {{-100, 100}, {100, 60}}, lineColor = {28, 108, 200}, horizontalAlignment = TextAlignment.Center, textString = "System Data")}),
|
|---|
| 304 | defaultComponentName = "SysData",
|
|---|
| 305 | defaultComponentPrefixes = "inner",
|
|---|
| 306 | missingInnerMessage = "
|
|---|
| 307 | No 'System Data' component is defined. A default component will be used, and generate a system base of 100 MVA, and a frequency of 50 Hz",
|
|---|
| 308 | Diagram(coordinateSystem(extent = {{-120, -100}, {120, 100}}, preserveAspectRatio = false)),
|
|---|
| 309 | defaultComponentPrefixes = "inner");
|
|---|
| 310 | end SystemBase;
|
|---|
| 311 |
|
|---|
| 312 | package BaseClasses
|
|---|
| 313 |
|
|---|
| 314 | partial model OnePort
|
|---|
| 315 | outer BugWithHomotopy.Elements.SystemBase SysData;
|
|---|
| 316 | parameter Types.ApparentPower SBase = SysData.SBase "System base power" annotation (Evaluate = true);
|
|---|
| 317 | parameter SI.Frequency fBase=SysData.fBase "System frequency" annotation ( Evaluate = true,Dialog(group="Power flow data"));
|
|---|
| 318 | parameter Types.Voltage VBase "Base voltage of the element" annotation (Evaluate = true, Dialog(group="Power flow data"));
|
|---|
| 319 |
|
|---|
| 320 | parameter Types.ActivePower P_0 "Initial active power" annotation (Dialog(group="Power flow data"));
|
|---|
| 321 | parameter Types.ReactivePower Q_0 "Initial reactive power" annotation (Dialog(group="Power flow data"));
|
|---|
| 322 | parameter Types.Voltage v_0 "Initial voltage magnitude (pu)" annotation (Dialog(group="Power flow data"));
|
|---|
| 323 | parameter Types.Angle angle_0 "Initial voltage angle" annotation (Dialog(group="Power flow data"));
|
|---|
| 324 |
|
|---|
| 325 | parameter Types.ComplexVoltage U_0 = Complex(v_0*cos(angle_0), v_0*sin(angle_0));
|
|---|
| 326 | parameter Types.ComplexPower S_0 = Complex(P_0, Q_0);
|
|---|
| 327 | parameter Types.ComplexCurrent I_0 = CM.conj(S_0/U_0);
|
|---|
| 328 | BugWithHomotopy.Interfaces.terminal p(
|
|---|
| 329 | v(re(nominal = VBase, start=U_0.re), im(nominal = VBase, start=U_0.im)),
|
|---|
| 330 | i(re(nominal = SBase/VBase, start=I_0.re), im(nominal = SBase/VBase, start=I_0.im))) annotation(
|
|---|
| 331 | Placement(visible = true, transformation(origin = {100, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {0, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
|
|---|
| 332 |
|
|---|
| 333 | Types.ActivePower P = S.re;
|
|---|
| 334 | Types.ReactivePower Q = S.im;
|
|---|
| 335 | Types.ComplexPower S=p.v*CM.conj(p.i);
|
|---|
| 336 | end OnePort;
|
|---|
| 337 |
|
|---|
| 338 | extends Modelica.Icons.BasesPackage;
|
|---|
| 339 |
|
|---|
| 340 | end BaseClasses;
|
|---|
| 341 | extends Modelica.Icons.Package;
|
|---|
| 342 |
|
|---|
| 343 | end Elements;
|
|---|
| 344 |
|
|---|
| 345 | package Types
|
|---|
| 346 |
|
|---|
| 347 | operator record ComplexPerUnit = Complex(re(unit = "1"), im(unit = "1"));
|
|---|
| 348 |
|
|---|
| 349 |
|
|---|
| 350 | operator record ComplexVoltage = SI.ComplexVoltage(re(nominal = 1e4, displayUnit="kV"), im(nominal = 1e4, displayUnit="kV"));
|
|---|
| 351 |
|
|---|
| 352 |
|
|---|
| 353 | operator record ComplexCurrent = SI.ComplexCurrent(re(nominal = 1e4, displayUnit="kA"), im(nominal = 1e4, displayUnit="kA"));
|
|---|
| 354 |
|
|---|
| 355 |
|
|---|
| 356 | operator record ComplexPower = SI.ComplexPower(re(nominal = 1e8, displayUnit="MW"), im(nominal = 1e8, displayUnit="MVA"));
|
|---|
| 357 |
|
|---|
| 358 |
|
|---|
| 359 | type ReactivePower = SI.ReactivePower(nominal = 1e6, displayUnit = "MVA");
|
|---|
| 360 |
|
|---|
| 361 |
|
|---|
| 362 | type ActivePower = SI.ActivePower(nominal = 1e6, displayUnit = "MW");
|
|---|
| 363 |
|
|---|
| 364 |
|
|---|
| 365 | type ApparentPower = SI.ApparentPower(nominal = 1e6, displayUnit = "MVA");
|
|---|
| 366 |
|
|---|
| 367 |
|
|---|
| 368 | type Angle = SI.Angle;
|
|---|
| 369 |
|
|---|
| 370 |
|
|---|
| 371 | type Voltage = SI.Voltage(nominal = 20e3, displayUnit = "kV");
|
|---|
| 372 |
|
|---|
| 373 | extends Modelica.Icons.TypesPackage;
|
|---|
| 374 |
|
|---|
| 375 | end Types;
|
|---|
| 376 |
|
|---|
| 377 | package Interfaces
|
|---|
| 378 |
|
|---|
| 379 | connector terminal "Connector for electrical blocks treating voltage and current as complex variables"
|
|---|
| 380 | Types.ComplexVoltage v(re(nominal=20e3, start=20e3), im(nominal=20e3, start=0)) "Voltage phasor";
|
|---|
| 381 | flow Types.ComplexCurrent i(re(nominal=1e6/20e3),im(nominal=1e6/20e3)) "Current phasor";
|
|---|
| 382 | annotation(
|
|---|
| 383 | Icon(graphics = {Rectangle(lineColor = {0, 85, 0}, fillColor = {0, 85, 0}, fillPattern = FillPattern.Solid, extent = {{-100, 100}, {100, -100}})}),
|
|---|
| 384 | Diagram(graphics = {Text(lineColor = {0, 0, 255}, extent = {{-100, 160}, {100, 120}}, textString = "%name"), Rectangle(lineColor = {0, 85, 0}, fillColor = {0, 85, 0}, fillPattern = FillPattern.Solid, extent = {{-100, 100}, {100, -100}})}));
|
|---|
| 385 | end terminal;
|
|---|
| 386 |
|
|---|
| 387 | extends Modelica.Icons.InterfacesPackage;
|
|---|
| 388 |
|
|---|
| 389 | end Interfaces;
|
|---|
| 390 | end BugWithHomotopy;
|
|---|