00001 #include "ElectricalModelEqns.h"
00002 using namespace std;
00003 using namespace adevs;
00004
00005 static const double pi = 3.1415926535897931;
00006
00007 #define TIME 0
00008 #define STATE_VARS_PER_GENR 5
00009 #define OMEGA(i) STATE_VARS_PER_GENR*i+1 // Angular speed
00010 #define THETA(i) STATE_VARS_PER_GENR*i+2 // Power angle
00011 #define C(i) STATE_VARS_PER_GENR*i+3 // Freqeuncy (mech. power) control signal
00012 #define Pm(i) STATE_VARS_PER_GENR*i+4 // Mechanical power
00013 #define EF(i) STATE_VARS_PER_GENR*i+5 // Exciter voltage
00014
00015 #define ALG_VARS_PER_BUS 2
00016 #define LoadRealCurrent(i) ALG_VARS_PER_BUS*i
00017 #define LoadImagCurrent(i) ALG_VARS_PER_BUS*i+1
00018
00019 #define EVENT_COND_PER_GENR 4
00020 #define OVER_FREQ_EVENT(i) (EVENT_COND_PER_GENR*i)
00021 #define UNDER_FREQ_EVENT(i) (EVENT_COND_PER_GENR*i+1)
00022 #define VOLTAGE_SATURATE_EVENT(i) (EVENT_COND_PER_GENR*i+2)
00023 #define VOLTAGE_UNSATURATE_EVENT(i) (EVENT_COND_PER_GENR*i+3)
00024 #define SENSOR_EVENT_START(num_genrs) (EVENT_COND_PER_GENR*num_genrs)
00025
00026
00027 static const double StateEventEps = 1E-5;
00028
00029 const int ElectricalModelEqns::SetLoad = -2;
00030
00031 const int ElectricalModelEqns::GenrTrip = -4;
00032
00033 const int ElectricalModelEqns::LoadAdj = -5;
00034
00035 const int ElectricalModelEqns::SensorOutput = -1;
00036
00037 const int ElectricalModelEqns::ACE = -3;
00038
00039 const int ElectricalModelEqns::LineAction = -6;
00040
00041
00042 ElectricalModelEqns::ElectricalModelEqns(ElectricalData* data, list<Sensor*>& sensors, bool no_events):
00043 dae_se1_system<PortValue<BasicEvent*> >(
00044
00045 data->getGenrCount()*STATE_VARS_PER_GENR+1,
00046
00047 (EVENT_COND_PER_GENR*data->getGenrCount()+sensors.size())*(!no_events),
00048
00049 data->getNodeCount()*ALG_VARS_PER_BUS
00050 ),
00051 num_nodes(data->getNodeCount()),
00052 Y(num_nodes)
00053 {
00054 init_model(data,sensors,no_events);
00055 }
00056
00057 ElectricalModelEqns::ElectricalModelEqns(ElectricalData* data, bool no_events):
00058 dae_se1_system<PortValue<BasicEvent*> >(
00059
00060 data->getGenrCount()*STATE_VARS_PER_GENR+1,
00061
00062 (EVENT_COND_PER_GENR*data->getGenrCount())*(!no_events),
00063
00064 data->getNodeCount()*ALG_VARS_PER_BUS
00065 ),
00066 num_nodes(data->getNodeCount()),
00067 Y(num_nodes)
00068 {
00069 list<Sensor*> empty_list;
00070 init_model(data,empty_list,no_events);
00071 }
00072
00073 void ElectricalModelEqns::init_model(ElectricalData* data, list<Sensor*>& sensors, bool no_events)
00074 {
00075 evar_iface = new SensorInterface(this);
00076
00077 this->data = data;
00078 this->no_events = no_events;
00079 this->load_adj = 0.0;
00080 Ef = new Complex[data->getGenrCount()];
00081
00082 breaker_closed = new bool[data->getGenrCount()];
00083
00084 excite_status = new ExcitationState[data->getGenrCount()];
00085
00086 data->buildAdmitMatrix(Y);
00087
00088 load_Yt = new Complex[num_nodes];
00089 genr_Yt = new Complex[num_nodes];
00090
00091 voltage = new Complex[num_nodes];
00092 current = new Complex[num_nodes];
00093 dV = new Complex[num_nodes];
00094 dI = new Complex[num_nodes];
00095 v0 = new Complex[num_nodes];
00096 load_freq = new double[num_nodes];
00097
00098 Pref = new double[data->getGenrCount()];
00099
00100 for (vector<unsigned>::const_iterator iter = data->getGenrs().begin();
00101 iter != data->getGenrs().end(); iter++)
00102 {
00103 genr_Yt[*iter] = 1.0/data->getGenrParams(*iter).Xd;
00104 }
00105 for (unsigned i = 0; i < num_nodes; i++)
00106 {
00107
00108 load_Yt[i] = data->getAdmittance(i);
00109 }
00110
00111 unsigned i;
00112 vector<unsigned>::const_iterator iter;
00113 for (i = 0, iter = data->getGenrs().begin();
00114 iter != data->getGenrs().end(); i++, iter++)
00115 {
00116
00117 genr_bus_lookup[*iter] = i;
00118
00119 ElectricalData::genr_t param = data->getGenrParams(*iter);
00120
00121 breaker_closed[i] = true;
00122
00123 excite_status[i] = NOT_SAT;
00124
00125 assert(fabs(param.Ef0) <= param.Ef_max);
00126 }
00127
00128 list<Sensor*>::iterator sensor_iter = sensors.begin();
00129 for (; sensor_iter != sensors.end(); sensor_iter++)
00130 {
00131 this->sensors.push_back(*sensor_iter);
00132 }
00133 }
00134
00135 void ElectricalModelEqns::init(double *q, double* a)
00136 {
00137
00138 q[TIME] = 0.0;
00139
00140 unsigned i;
00141 vector<unsigned>::const_iterator iter;
00142 for (i = 0, iter = data->getGenrs().begin();
00143 iter != data->getGenrs().end(); i++, iter++)
00144 {
00145
00146 ElectricalData::genr_t param = data->getGenrParams(*iter);
00147
00148 q[OMEGA(i)] = param.w0;
00149
00150 q[EF(i)] = param.Ef0;
00151
00152 q[Pm(i)] = param.Pm0;
00153 Pref[i] = param.Pe0;
00154 if (param.Pm0 != param.Pe0)
00155 {
00156 cerr << "WARNING: Genr " << i << " : Pm0-Pe0 = " << param.Pm0-param.Pe0 << endl;
00157 }
00158
00159 q[C(i)] = param.C0;
00160
00161 q[THETA(i)] = param.T0;
00162 }
00163
00164 for (i = 0; i < num_nodes; i++)
00165 current[i] = Complex(0.0,0.0);
00166
00167 for (i = 0, iter = data->getGenrs().begin();
00168 iter != data->getGenrs().end(); i++, iter++)
00169 {
00170
00171 Ef[i] = polar(q[EF(i)],q[THETA(i)]);
00172 current[*iter] += Ef[i]/data->getGenrParams(*iter).Xd;
00173 }
00174
00175 Y.solve_for_voltage(current,v0);
00176
00177 for (i = 0; i < num_nodes; i++)
00178 {
00179 load_freq[i] = 0.0;
00180 a[LoadRealCurrent(i)] = 0.0;
00181 a[LoadImagCurrent(i)] = 0.0;
00182 }
00183 }
00184
00185 void ElectricalModelEqns::solve_for_voltages(const double* q, const double* a)
00186 {
00187
00188 unsigned int i;
00189 vector<unsigned>::const_iterator iter;
00190
00191 for (i = 0; i < num_nodes; i++)
00192 current[i] = Complex(a[LoadRealCurrent(i)],a[LoadImagCurrent(i)]);
00193
00194 for (i = 0, iter = data->getGenrs().begin();
00195 iter != data->getGenrs().end(); i++, iter++)
00196 {
00197
00198 Ef[i] = polar(q[EF(i)],q[THETA(i)]);
00199
00200
00201 if (breaker_closed[i])
00202 current[*iter] += Ef[i]/data->getGenrParams(*iter).Xd;
00203 }
00204
00205 Y.solve_for_voltage(current,voltage);
00206 }
00207
00208 double ElectricalModelEqns::getVoltageWaveform(const double* q, unsigned bus)
00209 {
00216 double A[2][2], B[2], angle, mag, det;
00217
00218 Complex _Vgenr_;
00219 double _Fgenr_;
00220
00221 Complex _V_[num_nodes], _dV_[num_nodes];
00222 Complex _Iinj_[num_nodes];
00223
00224 double _Vt_ = 0.0;
00225 Complex VVV = 0.0;
00226 unsigned int i;
00227
00228 for (i = 0; i < num_nodes; i++)
00229 _Iinj_[i] = 0.0;
00230
00231 vector<unsigned>::const_iterator iter;
00232 for (i = 0, iter = data->getGenrs().begin();
00233 iter != data->getGenrs().end(); i++, iter++)
00234 {
00235
00236 if (breaker_closed[i])
00237 {
00238
00239 _Iinj_[*iter] = polar(q[EF(i)],q[THETA(i)])/data->getGenrParams(*iter).Xd;
00240 Y.solve_for_voltage(_Iinj_,_V_);
00241 _Vgenr_ = _V_[bus];
00242
00243 _Iinj_[*iter] = -polar(q[EF(i)]*q[OMEGA(i)],q[THETA(i)]-pi/2.0)/data->getGenrParams(*iter).Xd;
00244 Y.solve_for_voltage(_Iinj_,_dV_);
00245 angle = arg(_V_[bus]);
00246 mag = abs(_V_[bus]);
00247 A[0][0] = cos(angle);
00248 A[1][0] = -mag*cos(angle-pi/2.0);
00249 A[0][1] = sin(angle);
00250 A[1][1] = -mag*sin(angle-pi/2.0);
00251 B[0] = abs(_dV_[i])*cos(arg(_dV_[bus]));
00252 B[1] = abs(_dV_[i])*sin(arg(_dV_[bus]));
00253 det = A[0][0]*A[1][1]-A[1][0]*A[0][1];
00254 _Fgenr_ = (-A[0][1]*B[0]+A[0][0]*B[1])/det;
00255
00256 _Iinj_[*iter] = 0.0;
00257 }
00258 else
00259 {
00260 _Fgenr_ = 0.0;
00261 _Vgenr_ = 0.0;
00262 }
00263 _Vt_ += abs(_Vgenr_)*sin(2.0*pi*60.0*(1.0+_Fgenr_)*q[TIME]+arg(_Vgenr_));
00264
00265
00266 VVV += _Vgenr_;
00267 }
00268
00269 assert(fabs(real(VVV-voltage[bus])) < 1E-4 && fabs(imag(VVV-voltage[bus])) < 1E-4);
00270 return _Vt_;
00271 }
00272
00273 void ElectricalModelEqns::solve_for_freq(const double* q)
00274 {
00275
00276
00277
00278
00284 unsigned i;
00285 vector<unsigned>::const_iterator iter;
00286 for (i = 0; i < num_nodes; i++)
00287 dI[i] = Complex(0.0,0.0);
00288 for (i = 0, iter = data->getGenrs().begin();
00289 iter != data->getGenrs().end(); i++, iter++)
00290 {
00291 if (breaker_closed[i])
00292 {
00293 Complex V2 = polar(q[EF(i)]*q[OMEGA(i)],q[THETA(i)]-pi/2.0);
00294
00295
00296 dI[*iter] = -V2/data->getGenrParams(*iter).Xd;
00297 }
00298 }
00299 Y.solve_for_voltage(dI,dV);
00300 #pragma omp parallel for
00301 for (i = 0; i < num_nodes; i++)
00302 {
00303 double angle = arg(voltage[i]);
00304 double mag = abs(voltage[i]);
00305 double A[2][2], B[2];
00306 A[0][0] = cos(angle);
00307 A[1][0] = -mag*cos(angle-pi/2.0);
00308 A[0][1] = sin(angle);
00309 A[1][1] = -mag*sin(angle-pi/2.0);
00310 B[0] = abs(dV[i])*cos(arg(dV[i]));
00311 B[1] = abs(dV[i])*sin(arg(dV[i]));
00312 double det = A[0][0]*A[1][1]-A[1][0]*A[0][1];
00313 load_freq[i] = (-A[0][1]*B[0]+A[0][0]*B[1])/det;
00314 }
00315 }
00316
00317 void ElectricalModelEqns::alg_func(const double* q, const double* a, double* af)
00318 {
00319
00320 solve_for_voltages(q,a);
00321 solve_for_freq(q);
00329 unsigned i;
00330 #pragma omp parallel for
00331 for (i = 0; i < num_nodes; i++)
00332 {
00333 ElectricalData::load_param_t p(data->getLoadParams(i));
00334 double Kpz = 1.0-p.Kpi-p.Kpc-p.Kp1-p.Kp2;
00335 double Kqz = 1.0-p.Kqi-p.Kqc-p.Kq1-p.Kq2;
00336
00337
00338 if (abs(voltage[i]) <= 0.0)
00339 {
00340 af[LoadRealCurrent(i)] = 0.0;
00341 af[LoadImagCurrent(i)] = 0.0;
00342 }
00343 else
00344 {
00345 double v_v0 = abs(voltage[i])/abs(v0[i]);
00346 double freq = load_freq[i];
00347 double a_real =
00348 Kpz*v_v0*v_v0+
00349 p.Kpi*v_v0+
00350 p.Kpc+
00351 p.Kp1*pow(v_v0,p.npv1)*(1.0+p.npf1*freq)+
00352 p.Kp2*pow(v_v0,p.npv2)*(1.0+p.npf2*freq);
00353 double a_imag =
00354 Kqz*v_v0*v_v0+
00355 p.Kqi*v_v0+
00356 p.Kqc+
00357 p.Kq1*pow(v_v0,p.nqv1)*(1.0+p.nqf1*freq)+
00358 p.Kq2*pow(v_v0,p.nqv2)*(1.0+p.nqf2*freq);
00359
00360 Complex S0 = v0[i]*conj(v0[i]*load_Yt[i]*(1.0+load_adj));
00361 Complex S(a_real*real(S0),a_imag*imag(S0));
00362
00363 Complex Iinj = voltage[i]*load_Yt[i]*(1.0+load_adj) - conj(S/voltage[i]);
00364 af[LoadRealCurrent(i)] = real(Iinj);
00365 af[LoadImagCurrent(i)] = imag(Iinj);
00366 assert (!(af[LoadRealCurrent(i)] != af[LoadRealCurrent(i)]));
00367 assert (!(af[LoadImagCurrent(i)] != af[LoadImagCurrent(i)]));
00368 }
00369 }
00370 }
00371
00372 double ElectricalModelEqns::getBusFreq(unsigned bus)
00373 {
00374 return load_freq[bus];
00375 }
00376
00377 void ElectricalModelEqns::der_func(const double *q, const double* a, double *dq)
00378 {
00379
00380 solve_for_voltages(q,a);
00381
00382 unsigned i;
00383 vector<unsigned>::const_iterator iter;
00384 for (i = 0, iter = data->getGenrs().begin();
00385 iter != data->getGenrs().end(); i++, iter++)
00386 {
00387 Complex Vterm = voltage[*iter];
00388 ElectricalData::genr_t g = data->getGenrParams(*iter);
00389
00390 Complex I = (Ef[i]-Vterm)/g.Xd;
00391
00392 double Pe = real(Vterm*conj(I));
00393
00394 dq[OMEGA(i)] = ((q[Pm(i)]-Pe)/g.M-q[OMEGA(i)]*g.D)*
00395 breaker_closed[i];
00396 dq[THETA(i)] = q[OMEGA(i)]*breaker_closed[i];
00397
00398 dq[C(i)] = g.Tspd1*(Pref[i]-g.R*q[OMEGA(i)]-q[C(i)])*
00399 breaker_closed[i];
00400 dq[Pm(i)] = (g.Tspd2*(q[C(i)]-q[Pm(i)]))*breaker_closed[i];
00401
00402
00403
00404
00405 double Vt = abs(Vterm);
00406
00407 if (excite_status[i] == NOT_SAT || excite_status[i] == FALLING)
00408 dq[EF(i)] = ((g.Vref-Vt)*g.Te)*breaker_closed[i];
00409 else dq[EF(i)] = 0.0;
00410 }
00411
00412 dq[TIME] = 1.0;
00413 }
00414
00415 void ElectricalModelEqns::output_func(const double* q, const double* a, const bool* state_event,
00416 Bag<PortValue<BasicEvent*> >& yb)
00417 {
00418
00419 for (unsigned i = 0; i < sensors_primed.size(); i++)
00420 {
00421 yb.insert(PortValue<BasicEvent*>(SensorOutput,
00422 new SensorEvent(sensors_primed[i])));
00423 }
00424 }
00425
00426 void ElectricalModelEqns::internal_event(double* q, double* a, const bool* z)
00427 {
00428
00429 if (no_events) return;
00430
00431 unsigned i;
00432 vector<unsigned>::const_iterator iter;
00433 for (i = 0, iter = data->getGenrs().begin();
00434 iter != data->getGenrs().end(); i++, iter++)
00435 {
00436
00437 if (z[OVER_FREQ_EVENT(i)] || z[UNDER_FREQ_EVENT(i)])
00438 {
00439
00440 breaker_closed[i] = false;
00441
00442 Y.remove_self(*iter,genr_Yt[*iter]);
00443 genr_Yt[*iter] = Complex(0.0,0.0);
00444 }
00445
00446 if (z[VOLTAGE_SATURATE_EVENT(i)])
00447 {
00448 excite_status[i] = SAT;
00449 }
00450 else if (z[VOLTAGE_UNSATURATE_EVENT(i)])
00451 {
00452 if (excite_status[i] == SAT) excite_status[i] = FALLING;
00453 else excite_status[i] = NOT_SAT;
00454 }
00455 }
00456
00457 sensors_primed.clear();
00458
00459 evar_iface->setStateVars(q);
00460 for (i = 0; i < sensors.size(); i++)
00461 {
00462 int index = i + SENSOR_EVENT_START(data->getGenrCount());
00463 if (z[index] || (sensors[i]->getType() == Sensor::TIME_TR &&
00464 sensors[i]->distFromTrip(evar_iface) <= StateEventEps))
00465 {
00466 sensors_primed.push_back(sensors[i]);
00467 sensors[i]->trip(evar_iface);
00468 }
00469 }
00470 }
00471
00472 void ElectricalModelEqns::external_event(double* q, double* a, double e,
00473 const Bag<PortValue<BasicEvent*> >& xb)
00474 {
00475
00476 if (no_events) return;
00477
00478 for (Bag<PortValue<BasicEvent*> >::const_iterator iter = xb.begin();
00479 iter != xb.end(); iter++)
00480 {
00481
00482 if ((*iter).port == LoadAdj)
00483 {
00484 LoadAdjustEvent* load_event =
00485 dynamic_cast<LoadAdjustEvent*>((*iter).value);
00486 assert(load_event != NULL);
00487 for (unsigned bus = 0; bus < num_nodes; bus++)
00488 {
00489
00490 Y.remove_self(bus,(1.0+load_adj)*load_Yt[bus]);
00491
00492 Y.add_self(bus,(1.0+load_event->getAdjustment())*load_Yt[bus]);
00493 }
00494 load_adj = load_event->getAdjustment();
00495 }
00496
00497 if ((*iter).port == SetLoad)
00498 {
00499 LoadEvent* load_event = dynamic_cast<LoadEvent*>((*iter).value);
00500 assert(load_event != NULL);
00501 unsigned bus = load_event->getBusID();
00502
00503 Y.remove_self(bus,(1.0+load_adj)*load_Yt[bus]);
00504
00505 load_Yt[bus] = load_event->getAdmittance();
00506 Y.add_self(bus,(1.0+load_adj)*load_Yt[bus]);
00507 }
00508
00509 else if ((*iter).port == GenrTrip)
00510 {
00511 GenrFailEvent* genr_fail =
00512 dynamic_cast<GenrFailEvent*>((*iter).value);
00513 unsigned genr_index = genr_bus_lookup[genr_fail->getBusID()];
00514 breaker_closed[genr_index] = false;
00515
00516 Y.remove_self(genr_fail->getBusID(),genr_Yt[genr_fail->getBusID()]);
00517 genr_Yt[genr_fail->getBusID()] = Complex(0.0,0.0);
00518 }
00519
00520 else if ((*iter).port == LineAction)
00521 {
00522 LineEvent* line_event =
00523 dynamic_cast<LineEvent*>((*iter).value);
00524 ElectricalData::line_t which = line_event->getLine();
00525
00526 if (line_event->openLine())
00527 {
00528 if (which.tside >= 0)
00529 Y.remove_line(which.from,which.to,which.y,
00530 which.turns,which.phase_shift,which.tside);
00531 else Y.remove_line(which.from,which.to,which.y);
00532 }
00533
00534 else
00535 {
00536 if (which.tside >= 0)
00537 Y.add_line(which.from,which.to,which.y,
00538 which.turns,which.phase_shift,which.tside);
00539 else Y.add_line(which.from,which.to,which.y);
00540 }
00541 }
00542
00543 else if ((*iter).port == ACE)
00544 {
00545 unsigned i;
00546 vector<unsigned>::const_iterator iter;
00547 for (i = 0, iter = data->getGenrs().begin();
00548 iter != data->getGenrs().end(); i++, iter++)
00549 {
00550 ElectricalData::genr_t g = data->getGenrParams(*iter);
00551 Pref[i] -= g.Agc*q[OMEGA(i)];
00552 }
00553 }
00554 }
00555 }
00556
00557 void ElectricalModelEqns::confluent_event(double* q, double* a, const bool* z,
00558 const Bag<PortValue<BasicEvent*> >& xb)
00559 {
00560 internal_event(q,a,z);
00561 external_event(q,a,0,xb);
00562 }
00563
00564 void ElectricalModelEqns::state_event_func(const double* q, const double* a, double* z)
00565 {
00566
00567 if (no_events) return;
00568
00569 unsigned i;
00570 vector<unsigned>::const_iterator iter;
00571
00572 for (i = 0, iter = data->getGenrs().begin();
00573 iter != data->getGenrs().end(); i++, iter++)
00574 {
00575
00576 if (breaker_closed[i])
00577 {
00578
00579 z[OVER_FREQ_EVENT(i)] = 60.0*q[OMEGA(i)] -
00580 data->getGenrParams(*iter).FreqTol;
00581
00582 z[UNDER_FREQ_EVENT(i)] = data->getGenrParams(*iter).FreqTol +
00583 60.0*q[OMEGA(i)];
00584 }
00585 else z[OVER_FREQ_EVENT(i)] = z[UNDER_FREQ_EVENT(i)] = 1.0;
00586
00587
00588 if (excite_status[i] == NOT_SAT)
00589 {
00590 z[VOLTAGE_SATURATE_EVENT(i)] = q[EF(i)] -
00591 (data->getGenrParams(*iter).Ef_max+StateEventEps);
00592 z[VOLTAGE_UNSATURATE_EVENT(i)] = 1.0;
00593 }
00594
00595 else if (excite_status[i] == SAT)
00596 {
00597 z[VOLTAGE_SATURATE_EVENT(i)] = 1.0;
00598 z[VOLTAGE_UNSATURATE_EVENT(i)] = abs(voltage[*iter]) -
00599 (1.0+StateEventEps);
00600 }
00601 else if (excite_status[i] == FALLING)
00602 {
00603 z[VOLTAGE_SATURATE_EVENT(i)] = abs(voltage[*iter]) - 1.0;
00604 z[VOLTAGE_UNSATURATE_EVENT(i)] = q[EF(i)] -
00605 data->getGenrParams(*iter).Ef_max;
00606 }
00607 }
00608
00609 evar_iface->setStateVars(q);
00610 for (i = 0; i < sensors.size(); i++)
00611 {
00612 int index = i + SENSOR_EVENT_START(data->getGenrCount());
00613 if (sensors[i]->getType() == Sensor::STATE_TR)
00614 z[index] = sensors[i]->distFromTrip(evar_iface);
00615 else
00616 z[index] = -1.0;
00617 }
00618 }
00619
00620 double ElectricalModelEqns::time_event_func(const double* q, const double* a)
00621 {
00622
00623 if (!sensors_primed.empty()) return 0.0;
00624
00625 double h = DBL_MAX;
00626 if (!sensors.empty())
00627 {
00628 evar_iface->setStateVars(q);
00629 for (unsigned i = 0; i < sensors.size(); i++)
00630 {
00631 if (sensors[i]->getType() == Sensor::TIME_TR)
00632 {
00633 double htmp = sensors[i]->distFromTrip(evar_iface);
00634 if (htmp < h) h = htmp;
00635 }
00636 }
00637 }
00638 return h;
00639 }
00640
00641 void ElectricalModelEqns::gc_output(Bag<PortValue<BasicEvent*> >& gb)
00642 {
00643 for (Bag<PortValue<BasicEvent*> >::const_iterator iter = gb.begin();
00644 iter != gb.end(); iter++)
00645 {
00646 delete (*iter).value;
00647 }
00648 }
00649
00650 Complex ElectricalModelEqns::getVoltage(unsigned bus)
00651 {
00652 return voltage[bus];
00653 }
00654
00655 double ElectricalModelEqns::getGenrFreq(unsigned genr_number, const double *q)
00656 {
00657 return q[OMEGA(genr_number)];
00658 }
00659
00660 double ElectricalModelEqns::getRealGenrPower(unsigned genr_number, const double *q)
00661 {
00662 return q[Pm(genr_number)];
00663 }
00664
00665 bool ElectricalModelEqns::genrOffLine(unsigned genr_number)
00666 {
00667 return !breaker_closed[genr_number];
00668 }
00669
00670 Complex ElectricalModelEqns::getLoadPower(unsigned bus)
00671 {
00672 Complex V = getVoltage(bus);
00673 Complex Iinj(getAlgVar(LoadRealCurrent(bus)),
00674 getAlgVar(LoadImagCurrent(bus)));
00675 Complex y = load_Yt[bus]*(1.0+load_adj);
00676 Complex S = V*conj(V*y-Iinj);
00677 return S;
00678 }
00679
00680 void ElectricalModelEqns::postStep(double* q, double* a)
00681 {
00682
00683 solve_for_voltages(q,a);
00684 solve_for_freq(q);
00685 }
00686
00687 ElectricalModelEqns::~ElectricalModelEqns()
00688 {
00689 delete [] breaker_closed;
00690 delete [] voltage;
00691 delete [] current;
00692 delete [] dV;
00693 delete [] dI;
00694 delete [] v0;
00695 delete [] load_freq;
00696 delete [] load_Yt;
00697 delete [] genr_Yt;
00698 delete [] Ef;
00699 delete [] Pref;
00700 delete data;
00701 for (unsigned i = 0; i < sensors.size(); i++)
00702 delete sensors[i];
00703 }