00001
00002
00003
00004
00005
00006
00007
00008
00009
00012
00013
00015
00016
00017
00018
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "solver_nr.h"
00037
00038 #define MT // this enables multithreaded SuperLU
00039
00040 #ifdef MT
00041 #include <slu_mt_ddefs.h>
00042 #else
00043 #include <slu_ddefs.h>
00044 #endif
00045
00046
00047
00048 typedef struct {
00049 int *perm_c;
00050 int *perm_r;
00051 SuperMatrix A_LU;
00052 SuperMatrix B_LU;
00053 } SUPERLU_NR_vars;
00054
00055
00056 #include "powerflow.h"
00057
00058
00059 void sparse_init(SPARSE* sm, int nels, int ncols)
00060 {
00061 int indexval;
00062
00063
00064 sm->cols = (SP_E**)gl_malloc(ncols*sizeof(SP_E*));
00065
00066
00067 if (sm->cols == NULL)
00068 {
00069 GL_THROW("NR: Sparse matrix allocation failed");
00070
00071
00072
00073
00074 }
00075 else
00076 {
00077 for (indexval=0; indexval<ncols; indexval++)
00078 {
00079 sm->cols[indexval] = NULL;
00080 }
00081 }
00082
00083
00084 sm->llheap = (SP_E*)gl_malloc(nels*sizeof(SP_E));
00085
00086
00087 if (sm->llheap == NULL)
00088 {
00089 GL_THROW("NR: Sparse matrix allocation failed");
00090
00091 }
00092 else
00093 {
00094 for (indexval=0; indexval<nels; indexval++)
00095 {
00096 sm->llheap[indexval].next = NULL;
00097 sm->llheap[indexval].row_ind = -1;
00098 sm->llheap[indexval].value = 0.0;
00099 }
00100 }
00101
00102
00103 sm->llptr = 0;
00104 sm->ncols = ncols;
00105 }
00106
00107
00108 void sparse_clear(SPARSE* sm)
00109 {
00110
00111 gl_free(sm->llheap);
00112 gl_free(sm->cols);
00113
00114
00115 sm->llheap = NULL;
00116 sm->cols = NULL;
00117
00118
00119 sm->ncols = 0;
00120 }
00121
00122 void sparse_reset(SPARSE* sm, int ncols)
00123 {
00124 int indexval;
00125
00126
00127 sm->ncols = ncols;
00128
00129
00130 for (indexval=0; indexval<ncols; indexval++)
00131 {
00132 sm->cols[indexval] = NULL;
00133 }
00134
00135
00136 sm->llptr = 0;
00137 }
00138
00139
00140 inline void sparse_add(SPARSE* sm, int row, int col, double value, BUSDATA *bus_values, unsigned int bus_values_count, NR_SOLVER_STRUCT *powerflow_information, int island_number_curr)
00141 {
00142 unsigned int bus_index_val, bus_start_val, bus_end_val;
00143 bool found_proper_bus_val;
00144
00145 SP_E* insertion_point = sm->cols[col];
00146 SP_E* new_list_element = &(sm->llheap[sm->llptr++]);
00147
00148 new_list_element->next = NULL;
00149 new_list_element->row_ind = row;
00150 new_list_element->value = value;
00151
00152
00153 if(insertion_point != NULL)
00154 {
00155 if(insertion_point->row_ind > new_list_element->row_ind)
00156 {
00157
00158 new_list_element->next = insertion_point;
00159 sm->cols[col] = new_list_element;
00160 }
00161 else
00162 {
00163 while((insertion_point->next != NULL) && (insertion_point->next->row_ind < new_list_element->row_ind))
00164 {
00165 insertion_point = insertion_point->next;
00166 }
00167
00168
00169 if (insertion_point->next != NULL)
00170 {
00171 if (insertion_point->next->row_ind == new_list_element->row_ind)
00172 {
00173
00174 found_proper_bus_val = false;
00175
00176
00177 for (bus_index_val=0; bus_index_val<bus_values_count; bus_index_val++)
00178 {
00179
00180 if (bus_values[bus_index_val].island_number == island_number_curr)
00181 {
00182
00183 bus_start_val = 2*bus_values[bus_index_val].Matrix_Loc;
00184 bus_end_val = bus_start_val + 2*powerflow_information->BA_diag[bus_index_val].size - 1;
00185
00186
00187 if ((new_list_element->row_ind >= bus_start_val) && (new_list_element->row_ind <= bus_end_val))
00188 {
00189
00190 if (bus_values[bus_index_val].name != NULL)
00191 {
00192 GL_THROW("NR: duplicate admittance entry found - attaches to node %s - check for parallel circuits between common nodes!",bus_values[bus_index_val].name);
00193
00194
00195
00196
00197
00198 }
00199 else
00200 {
00201 GL_THROW("NR: duplicate admittance entry found - no name available - check for parallel circuits between common nodes!");
00202
00203
00204
00205
00206
00207 }
00208 }
00209
00210 }
00211
00212 }
00213 }
00214 }
00215 else
00216 {
00217 if (insertion_point->row_ind == new_list_element->row_ind)
00218 {
00219
00220 found_proper_bus_val = false;
00221
00222
00223 for (bus_index_val=0; bus_index_val<bus_values_count; bus_index_val++)
00224 {
00225
00226 if (bus_values[bus_index_val].island_number == island_number_curr)
00227 {
00228
00229 bus_start_val = 2*bus_values[bus_index_val].Matrix_Loc;
00230 bus_end_val = bus_start_val + 2*powerflow_information->BA_diag[bus_index_val].size - 1;
00231
00232
00233 if ((new_list_element->row_ind >= bus_start_val) && (new_list_element->row_ind <= bus_end_val))
00234 {
00235
00236 if (bus_values[bus_index_val].name != NULL)
00237 {
00238 GL_THROW("NR: duplicate admittance entry found - attaches to node %s - check for parallel circuits between common nodes!",bus_values[bus_index_val].name);
00239
00240 }
00241 else
00242 {
00243 GL_THROW("NR: duplicate admittance entry found - no name available - check for parallel circuits between common nodes!");
00244
00245 }
00246 }
00247
00248 }
00249
00250 }
00251 }
00252 }
00253
00254
00255 new_list_element->next = insertion_point->next;
00256 insertion_point->next = new_list_element;
00257 }
00258 }
00259 else
00260 sm->cols[col] = new_list_element;
00261 }
00262
00263 void sparse_tonr(SPARSE* sm, NR_SOLVER_VARS *matrices_LUin)
00264 {
00265
00266 unsigned int rowidx = 0;
00267 unsigned int colidx = 0;
00268 unsigned int i;
00269 SP_E* LL_pointer;
00270
00271 matrices_LUin->cols_LU[0] = 0;
00272 LL_pointer = NULL;
00273 for(i = 0; i < sm->ncols; i++)
00274 {
00275 LL_pointer = sm->cols[i];
00276 if(LL_pointer != NULL)
00277 {
00278 matrices_LUin->cols_LU[colidx++] = rowidx;
00279 }
00280 while(LL_pointer != NULL)
00281 {
00282 matrices_LUin->rows_LU[rowidx] = LL_pointer->row_ind;
00283 matrices_LUin->a_LU[rowidx] = LL_pointer->value;
00284 ++rowidx;
00285 LL_pointer = LL_pointer->next;
00286 }
00287 }
00288 }
00289
00297 int64 solver_nr(unsigned int bus_count, BUSDATA *bus, unsigned int branch_count, BRANCHDATA *branch, NR_SOLVER_STRUCT *powerflow_values, NRSOLVERMODE powerflow_type , NR_MESHFAULT_IMPEDANCE *mesh_imped_vals, bool *bad_computations)
00298 {
00299
00300 int island_index_val;
00301
00302
00303 int island_loop_index;
00304
00305
00306 FILE *FPoutVal;
00307
00308
00309 int func_result_val;
00310
00311
00312 FUNCTIONADDR temp_fxn_val;
00313
00314
00315 STATUS call_return_status;
00316
00317
00318 unsigned char phase_worka, phase_workb, phase_workc, phase_workd, phase_worke;
00319
00320
00321 double tempIcalcReal, tempIcalcImag;
00322 double tempPbus;
00323 double tempQbus;
00324
00325
00326 unsigned int indexer, tempa, tempb, jindexer, kindexer;
00327 char jindex, kindex;
00328 char temp_index, temp_index_b;
00329 unsigned int temp_index_c;
00330
00331
00332 complex tempY[3][3];
00333
00334
00335 double temp_z_store[6][6];
00336
00337
00338 bool Full_Mat_A, Full_Mat_B, proceed_flag;
00339
00340
00341 complex temp_complex_0, temp_complex_1, temp_complex_2, temp_complex_3, temp_complex_4, temp_complex_5;
00342 complex aval, avalsq;
00343
00344
00345 char temp_size, temp_size_b, temp_size_c;
00346
00347
00348 complex Temp_Ad_A[3][3];
00349 complex Temp_Ad_B[3][3];
00350
00351
00352 complex DVConvCheck[3];
00353 double CurrConvVal;
00354
00355
00356 double work_vals_double_0, work_vals_double_1,work_vals_double_2,work_vals_double_3,work_vals_double_4;
00357 char work_vals_char_0;
00358 bool something_has_been_output;
00359
00360
00361 SuperMatrix L_LU,U_LU;
00362 NCformat *Astore;
00363 DNformat *Bstore;
00364 int nnz, info;
00365 unsigned int m,n;
00366 double *sol_LU;
00367
00368
00369 SP_E *temp_element;
00370 int row, col;
00371 double value;
00372
00373
00374 bool still_iterating_islands;
00375 bool proceed_to_next_island;
00376 int64 return_value_for_solver_NR;
00377
00378
00379 SUPERLU_NR_vars *curr_island_superLU_vars;
00380
00381 #ifndef MT
00382 superlu_options_t options;
00383 SuperLUStat_t stat;
00384 #endif
00385
00386
00387 NR_solver_working = true;
00388
00389
00390 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00391 {
00392 powerflow_values->island_matrix_values[island_loop_index].index_count = 0;
00393 }
00394
00395
00396 *bad_computations = false;
00397
00398
00399 if (powerflow_type != PF_NORMAL)
00400 {
00401 if (powerflow_type == PF_DYNCALC)
00402 {
00403
00404 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00405 {
00406 powerflow_values->island_matrix_values[island_loop_index].swing_is_a_swing = false;
00407 }
00408
00409
00410 for (indexer=0; indexer<bus_count; indexer++)
00411 {
00412
00413 if ((bus[indexer].type > 1) && (bus[indexer].swing_functions_enabled == true))
00414 {
00415
00416 if ((*bus[indexer].dynamics_enabled==true) && (bus[indexer].DynCurrent != NULL))
00417 {
00418
00419 bus[indexer].swing_functions_enabled = false;
00420 }
00421 }
00422
00423 }
00424 }
00425 else
00426 {
00427
00428 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00429 {
00430
00431 powerflow_values->island_matrix_values[island_loop_index].swing_is_a_swing = true;
00432 }
00433 }
00434 }
00435 else
00436 {
00437
00438 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00439 {
00440 powerflow_values->island_matrix_values[island_loop_index].swing_is_a_swing = true;
00441 }
00442 }
00443
00444
00445 if (powerflow_type == PF_DYNINIT)
00446 {
00447
00448 aval = complex(-0.5,(sqrt(3.0)/2.0));
00449 avalsq = aval*aval;
00450 }
00451 else
00452 {
00453 aval = 0.0;
00454 avalsq = 0.0;
00455 }
00456
00457 if (matrix_solver_method==MM_EXTERN)
00458 {
00459 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00460 {
00461
00462 powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars = ((void *(*)(void *))(LUSolverFcns.ext_init))(powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars);
00463
00464
00465 if (powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars==NULL)
00466 {
00467 GL_THROW("External LU matrix solver failed to allocate memory properly!");
00468
00469
00470
00471
00472
00473 }
00474 }
00475 }
00476 else if (matrix_solver_method == MM_SUPERLU)
00477 {
00478
00479 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00480 {
00481
00482 if (powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars == NULL)
00483 {
00484
00485 curr_island_superLU_vars = (SUPERLU_NR_vars *)gl_malloc(sizeof(SUPERLU_NR_vars));
00486
00487
00488 if (curr_island_superLU_vars == NULL)
00489 {
00490 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
00491
00492 }
00493
00494
00495 curr_island_superLU_vars->A_LU.Store = NULL;
00496
00497
00498 curr_island_superLU_vars->A_LU.Stype = SLU_NC;
00499 curr_island_superLU_vars->A_LU.Dtype = SLU_D;
00500 curr_island_superLU_vars->A_LU.Mtype = SLU_GE;
00501 curr_island_superLU_vars->A_LU.nrow = 0;
00502 curr_island_superLU_vars->A_LU.ncol = 0;
00503
00504
00505 curr_island_superLU_vars->B_LU.Store = NULL;
00506
00507
00508 curr_island_superLU_vars->B_LU.Stype = SLU_NC;
00509 curr_island_superLU_vars->B_LU.Dtype = SLU_D;
00510 curr_island_superLU_vars->B_LU.Mtype = SLU_GE;
00511 curr_island_superLU_vars->B_LU.nrow = 0;
00512 curr_island_superLU_vars->B_LU.ncol = 0;
00513
00514
00515 curr_island_superLU_vars->perm_c = NULL;
00516 curr_island_superLU_vars->perm_r = NULL;
00517
00518
00519 powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars = (void *)curr_island_superLU_vars;
00520 }
00521 }
00522
00523
00524 curr_island_superLU_vars = NULL;
00525 }
00526
00527 if (NR_admit_change)
00528 {
00529
00530 if (powerflow_values->BA_diag == NULL)
00531 {
00532 powerflow_values->BA_diag = (Bus_admit *)gl_malloc(bus_count *sizeof(Bus_admit));
00533
00534
00535 if (powerflow_values->BA_diag == NULL)
00536 {
00537 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
00538
00539
00540
00541
00542
00543 }
00544 }
00545
00546
00547 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00548 {
00549 powerflow_values->island_matrix_values[island_loop_index].bus_count = 0;
00550 }
00551
00552 for (indexer=0; indexer<bus_count; indexer++)
00553 {
00554
00555 if (bus[indexer].island_number != -1)
00556 {
00557 powerflow_values->island_matrix_values[bus[indexer].island_number].bus_count++;
00558 }
00559
00560
00561 if ((bus[indexer].phases & 0x80) == 0x80)
00562 powerflow_values->BA_diag[indexer].size = 2;
00563 else
00564 {
00565 phase_worka = 0;
00566 for (jindex=0; jindex<3; jindex++)
00567 {
00568 phase_worka += ((bus[indexer].phases & (0x01 << jindex)) >> jindex);
00569 }
00570 powerflow_values->BA_diag[indexer].size = phase_worka;
00571 }
00572
00573
00574 for (jindex=0; jindex<3; jindex++)
00575 {
00576 for (kindex=0; kindex<3; kindex++)
00577 {
00578 powerflow_values->BA_diag[indexer].Y[jindex][kindex] = 0;
00579 tempY[jindex][kindex] = 0;
00580 }
00581 }
00582
00583
00584 if ((powerflow_type!=PF_NORMAL) && (bus[indexer].full_Y != NULL))
00585 {
00586
00587 for (jindex=0; jindex<3; jindex++)
00588 {
00589 for (kindex=0; kindex<3; kindex++)
00590 {
00591 tempY[jindex][kindex] = bus[indexer].full_Y[jindex*3+kindex];
00592 }
00593 }
00594 }
00595
00596
00597 for (kindexer=0; kindexer<(bus[indexer].Link_Table_Size);kindexer++)
00598 {
00599
00600 jindexer = bus[indexer].Link_Table[kindexer];
00601
00602 if ((branch[jindexer].from == indexer) || (branch[jindexer].to == indexer))
00603 {
00604 if ((bus[indexer].phases & 0x07) == 0x07)
00605 {
00606 for (jindex=0; jindex<3; jindex++)
00607 {
00608
00609 phase_workb = 0x04 >> jindex;
00610
00611 if ((phase_workb & branch[jindexer].phases) == phase_workb)
00612 {
00613 for (kindex=0; kindex<3; kindex++)
00614 {
00615
00616 phase_workd = 0x04 >> kindex;
00617
00618 if ((phase_workd & branch[jindexer].phases) == phase_workd)
00619 {
00620 if (branch[jindexer].from == indexer)
00621 {
00622 tempY[jindex][kindex] += branch[jindexer].YSfrom[jindex*3+kindex];
00623 }
00624 else
00625 {
00626 tempY[jindex][kindex] += branch[jindexer].YSto[jindex*3+kindex];
00627 }
00628 }
00629 }
00630 }
00631 }
00632 }
00633 else if ((bus[indexer].phases & 0x80) == 0x80)
00634 {
00635 if (branch[jindexer].from == indexer)
00636 {
00637
00638 if ((bus[indexer].phases & 0x20) == 0x20)
00639 {
00640
00641 tempY[0][0] -= branch[jindexer].YSfrom[0];
00642 tempY[0][1] -= branch[jindexer].YSfrom[1];
00643 tempY[1][0] -= branch[jindexer].YSfrom[3];
00644 tempY[1][1] -= branch[jindexer].YSfrom[4];
00645 }
00646 else
00647 {
00648 tempY[0][0] += branch[jindexer].YSfrom[0];
00649 tempY[0][1] += branch[jindexer].YSfrom[1];
00650 tempY[1][0] += branch[jindexer].YSfrom[3];
00651 tempY[1][1] += branch[jindexer].YSfrom[4];
00652 }
00653 }
00654 else
00655 {
00656
00657
00658 if (((bus[indexer].phases & 0x20) == 0x20) && (branch[jindexer].lnk_type == 1))
00659 {
00660 tempY[0][0] -= branch[jindexer].YSto[0];
00661 tempY[0][1] -= branch[jindexer].YSto[1];
00662 tempY[1][0] -= branch[jindexer].YSto[3];
00663 tempY[1][1] -= branch[jindexer].YSto[4];
00664 }
00665 else
00666 {
00667 tempY[0][0] += branch[jindexer].YSto[0];
00668 tempY[0][1] += branch[jindexer].YSto[1];
00669 tempY[1][0] += branch[jindexer].YSto[3];
00670 tempY[1][1] += branch[jindexer].YSto[4];
00671 }
00672 }
00673 }
00674 else
00675 {
00676 switch(bus[indexer].phases & 0x07) {
00677 case 0x00:
00678 {
00679 break;
00680 }
00681 case 0x01:
00682 {
00683 if ((branch[jindexer].phases & 0x01) == 0x01)
00684 {
00685 if (branch[jindexer].from == indexer)
00686 {
00687 tempY[0][0] += branch[jindexer].YSfrom[8];
00688 }
00689 else
00690 {
00691 tempY[0][0] += branch[jindexer].YSto[8];
00692 }
00693 }
00694 break;
00695 }
00696 case 0x02:
00697 {
00698 if ((branch[jindexer].phases & 0x02) == 0x02)
00699 {
00700 if (branch[jindexer].from == indexer)
00701 {
00702 tempY[0][0] += branch[jindexer].YSfrom[4];
00703 }
00704 else
00705 {
00706 tempY[0][0] += branch[jindexer].YSto[4];
00707 }
00708 }
00709 break;
00710 }
00711 case 0x03:
00712 {
00713 phase_worka = (branch[jindexer].phases & 0x03);
00714
00715 if (phase_worka == 0x03)
00716 {
00717 if (branch[jindexer].from == indexer)
00718 {
00719 tempY[0][0] += branch[jindexer].YSfrom[4];
00720 tempY[0][1] += branch[jindexer].YSfrom[5];
00721 tempY[1][0] += branch[jindexer].YSfrom[7];
00722 tempY[1][1] += branch[jindexer].YSfrom[8];
00723 }
00724 else
00725 {
00726 tempY[0][0] += branch[jindexer].YSto[4];
00727 tempY[0][1] += branch[jindexer].YSto[5];
00728 tempY[1][0] += branch[jindexer].YSto[7];
00729 tempY[1][1] += branch[jindexer].YSto[8];
00730 }
00731 }
00732 else if (phase_worka == 0x01)
00733 {
00734 if (branch[jindexer].from == indexer)
00735 {
00736 tempY[1][1] += branch[jindexer].YSfrom[8];
00737 }
00738 else
00739 {
00740 tempY[1][1] += branch[jindexer].YSto[8];
00741 }
00742 }
00743 else if (phase_worka == 0x02)
00744 {
00745 if (branch[jindexer].from == indexer)
00746 {
00747 tempY[0][0] += branch[jindexer].YSfrom[4];
00748 }
00749 else
00750 {
00751 tempY[0][0] += branch[jindexer].YSto[4];
00752 }
00753 }
00754 else
00755 ;
00756 break;
00757 }
00758 case 0x04:
00759 {
00760 if ((branch[jindexer].phases & 0x04) == 0x04)
00761 {
00762 if (branch[jindexer].from == indexer)
00763 {
00764 tempY[0][0] += branch[jindexer].YSfrom[0];
00765 }
00766 else
00767 {
00768 tempY[0][0] += branch[jindexer].YSto[0];
00769 }
00770 }
00771 break;
00772 }
00773 case 0x05:
00774 {
00775 phase_worka = branch[jindexer].phases & 0x05;
00776
00777 if (phase_worka == 0x05)
00778 {
00779 if (branch[jindexer].from == indexer)
00780 {
00781 tempY[0][0] += branch[jindexer].YSfrom[0];
00782 tempY[0][1] += branch[jindexer].YSfrom[2];
00783 tempY[1][0] += branch[jindexer].YSfrom[6];
00784 tempY[1][1] += branch[jindexer].YSfrom[8];
00785 }
00786 else
00787 {
00788 tempY[0][0] += branch[jindexer].YSto[0];
00789 tempY[0][1] += branch[jindexer].YSto[2];
00790 tempY[1][0] += branch[jindexer].YSto[6];
00791 tempY[1][1] += branch[jindexer].YSto[8];
00792 }
00793 }
00794 else if (phase_worka == 0x04)
00795 {
00796 if (branch[jindexer].from == indexer)
00797 {
00798 tempY[0][0] += branch[jindexer].YSfrom[0];
00799 }
00800 else
00801 {
00802 tempY[0][0] += branch[jindexer].YSto[0];
00803 }
00804 }
00805 else if (phase_worka == 0x01)
00806 {
00807 if (branch[jindexer].from == indexer)
00808 {
00809 tempY[1][1] += branch[jindexer].YSfrom[8];
00810 }
00811 else
00812 {
00813 tempY[1][1] += branch[jindexer].YSto[8];
00814 }
00815 }
00816 else
00817 ;
00818 break;
00819 }
00820 case 0x06:
00821 {
00822 phase_worka = branch[jindexer].phases & 0x06;
00823
00824 if (phase_worka == 0x06)
00825 {
00826 if (branch[jindexer].from == indexer)
00827 {
00828 tempY[0][0] += branch[jindexer].YSfrom[0];
00829 tempY[0][1] += branch[jindexer].YSfrom[1];
00830 tempY[1][0] += branch[jindexer].YSfrom[3];
00831 tempY[1][1] += branch[jindexer].YSfrom[4];
00832 }
00833 else
00834 {
00835 tempY[0][0] += branch[jindexer].YSto[0];
00836 tempY[0][1] += branch[jindexer].YSto[1];
00837 tempY[1][0] += branch[jindexer].YSto[3];
00838 tempY[1][1] += branch[jindexer].YSto[4];
00839 }
00840 }
00841 else if (phase_worka == 0x04)
00842 {
00843 if (branch[jindexer].from == indexer)
00844 {
00845 tempY[0][0] += branch[jindexer].YSfrom[0];
00846 }
00847 else
00848 {
00849 tempY[0][0] += branch[jindexer].YSto[0];
00850 }
00851 }
00852 else if (phase_worka == 0x02)
00853 {
00854 if (branch[jindexer].from == indexer)
00855 {
00856 tempY[1][1] += branch[jindexer].YSfrom[4];
00857 }
00858 else
00859 {
00860 tempY[1][1] += branch[jindexer].YSto[4];
00861 }
00862 }
00863 else
00864 ;
00865 break;
00866 }
00867 default:
00868 {
00869 GL_THROW("Unknown phase connection in NR self admittance diagonal");
00870
00871
00872
00873
00874
00875 break;
00876 }
00877 }
00878 }
00879 }
00880 else
00881 ;
00882 }
00883
00884
00885
00886 island_index_val = bus[indexer].island_number;
00887
00888 if (island_index_val != -1)
00889 {
00890 powerflow_values->BA_diag[indexer].col_ind = powerflow_values->BA_diag[indexer].row_ind = powerflow_values->island_matrix_values[island_index_val].index_count;
00891 bus[indexer].Matrix_Loc = powerflow_values->island_matrix_values[island_index_val].index_count;
00892 powerflow_values->island_matrix_values[island_index_val].index_count += powerflow_values->BA_diag[indexer].size;
00893 }
00894
00895
00896 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
00897 {
00898 for (kindex=0; kindex<powerflow_values->BA_diag[indexer].size; kindex++)
00899 {
00900 powerflow_values->BA_diag[indexer].Y[jindex][kindex] = tempY[jindex][kindex];
00901 }
00902 }
00903
00904
00905 if (bus[indexer].full_Y_all != NULL)
00906 {
00907 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
00908 {
00909 for (kindex=0; kindex<powerflow_values->BA_diag[indexer].size; kindex++)
00910 {
00911 bus[indexer].full_Y_all[jindex*3+kindex] = tempY[jindex][kindex];
00912 }
00913 }
00914 }
00915 }
00916
00917
00918 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
00919 {
00920
00921 powerflow_values->island_matrix_values[island_loop_index].total_variables = powerflow_values->island_matrix_values[island_loop_index].index_count;
00922
00923
00924 if (powerflow_values->island_matrix_values[island_loop_index].total_variables > powerflow_values->island_matrix_values[island_loop_index].max_total_variables)
00925 powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed = true;
00926
00928
00929
00930 powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ = 0;
00931 }
00932
00933 for (jindexer=0; jindexer<branch_count;jindexer++)
00934 {
00935 tempa = branch[jindexer].from;
00936 tempb = branch[jindexer].to;
00937
00938
00939 island_index_val = bus[tempa].island_number;
00940
00941
00942 if (island_index_val == -1)
00943 {
00944 continue;
00945 }
00946
00947
00948 if ((bus[tempa].Matrix_Loc == -1) || (bus[tempb].Matrix_Loc == -1))
00949 {
00950 GL_THROW("An element in NR line:%d was not properly localized");
00951
00952
00953
00954
00955
00956 }
00957
00958 if (((branch[jindexer].phases & 0x80) == 0x80) && (branch[jindexer].v_ratio==1.0))
00959 {
00960 for (jindex=0; jindex<2; jindex++)
00961 {
00962 for (kindex=0; kindex<2; kindex++)
00963 {
00964 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
00965 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
00966
00967 if (((branch[jindexer].Yto[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
00968 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
00969
00970 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
00971 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
00972
00973 if (((branch[jindexer].Yto[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
00974 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
00975 }
00976 }
00977 }
00978 else
00979 {
00980
00981 if ((branch[jindexer].phases & 0x80) != 0x80)
00982 {
00983 for (jindex=0; jindex<3; jindex++)
00984 {
00985
00986 phase_workb = 0x04 >> jindex;
00987
00988 if ((phase_workb & branch[jindexer].phases) == phase_workb)
00989 {
00990 for (kindex=0; kindex<3; kindex++)
00991 {
00992
00993 phase_workd = 0x04 >> kindex;
00994
00995 if ((phase_workd & branch[jindexer].phases) == phase_workd)
00996 {
00997 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
00998 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
00999
01000 if (((branch[jindexer].Yto[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01001 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
01002
01003 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01004 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
01005
01006 if (((branch[jindexer].Yto[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01007 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
01008 }
01009 }
01010 }
01011 }
01012 }
01013 else
01014 {
01015 for (jindex=0; jindex<3; jindex++)
01016 {
01017
01018 phase_workb = 0x04 >> jindex;
01019
01020 if ((phase_workb & branch[jindexer].phases) == phase_workb)
01021 {
01022 for (kindex=0; kindex<3; kindex++)
01023 {
01024 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01025 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
01026
01027 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01028 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
01029 }
01030
01031
01032 for (kindex=0; kindex<3; kindex++)
01033 {
01034 if (((branch[jindexer].Yto[kindex*3+jindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01035 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
01036
01037 if (((branch[jindexer].Yto[kindex*3+jindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01038 powerflow_values->island_matrix_values[island_index_val].size_offdiag_PQ += 1;
01039 }
01040 }
01041 }
01042 }
01043 }
01044 }
01045
01046
01047 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
01048 {
01049
01050 if (powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ == NULL)
01051 {
01052 powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ = (Y_NR *)gl_malloc((powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ*2) *sizeof(Y_NR));
01053
01054
01055 if (powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ == NULL)
01056 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
01057
01058
01059 powerflow_values->island_matrix_values[island_loop_index].max_size_offdiag_PQ = powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ;
01060 }
01061 else if (powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ > powerflow_values->island_matrix_values[island_loop_index].max_size_offdiag_PQ)
01062 {
01063
01064 gl_free(powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ);
01065
01066
01067 powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ = (Y_NR *)gl_malloc((powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ*2) *sizeof(Y_NR));
01068
01069
01070 if (powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ == NULL)
01071 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
01072
01073
01074 powerflow_values->island_matrix_values[island_loop_index].max_size_offdiag_PQ = powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ;
01075
01076
01077 powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed = true;
01078 }
01079
01080
01081 powerflow_values->island_matrix_values[island_loop_index].indexer = 0;
01082 }
01083
01084 for (jindexer=0; jindexer<branch_count;jindexer++)
01085 {
01086
01087 tempa = branch[jindexer].from;
01088 tempb = branch[jindexer].to;
01089
01090
01091 island_index_val = bus[tempa].island_number;
01092
01093
01094 if (island_index_val == -1)
01095 {
01096 continue;
01097 }
01098
01099 phase_worka = 0;
01100 phase_workb = 0;
01101 for (jindex=0; jindex<3; jindex++)
01102 {
01103 phase_worka += ((bus[tempa].phases & (0x01 << jindex)) >> jindex);
01104 phase_workb += ((bus[tempb].phases & (0x01 << jindex)) >> jindex);
01105 }
01106
01107 if ((phase_worka==3) && (phase_workb==3))
01108 {
01109 for (jindex=0; jindex<3; jindex++)
01110 {
01111
01112 phase_workd = 0x04 >> jindex;
01113
01114 if ((branch[jindexer].phases & phase_workd) == phase_workd)
01115 {
01116 for (kindex=0; kindex<3; kindex++)
01117 {
01118
01119 phase_worke = 0x04 >> kindex;
01120
01121 if ((branch[jindexer].phases & phase_worke) == phase_worke)
01122 {
01123
01124 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01125 {
01126 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01127 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01128 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Im());
01129 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01130 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 3;
01131 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 3;
01132 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yfrom[jindex*3+kindex]).Im();
01133 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01134 }
01135
01136 if (((branch[jindexer].Yto[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01137 {
01138 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01139 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01140 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Im());
01141 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01142 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 3;
01143 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 3;
01144 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yto[jindex*3+kindex]).Im();
01145 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01146 }
01147
01148 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01149 {
01150 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 3;
01151 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01152 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01153 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01154 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01155 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 3;
01156 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01157 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01158 }
01159
01160 if (((branch[jindexer].Yto[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01161 {
01162 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 3;
01163 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01164 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Re());
01165 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01166 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01167 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 3;
01168 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Re());
01169 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01170 }
01171 }
01172 }
01173 }
01174 }
01175 }
01176 else if (((bus[tempa].phases & 0x80) == 0x80) || ((bus[tempb].phases & 0x80) == 0x80))
01177 {
01178 if (((bus[tempa].phases & 0x80) == 0x80) && ((bus[tempb].phases & 0x80) == 0x80))
01179 {
01180 for (jindex=0; jindex<2; jindex++)
01181 {
01182 for (kindex=0; kindex<2; kindex++)
01183 {
01184
01185 if (((bus[tempa].phases & 0x20) & (bus[tempb].phases & 0x20)) == 0x20)
01186 {
01187 GL_THROW("NR: SPCT to SPCT via triplex connections are unsupported at this time.");
01188
01189
01190
01191
01192
01193 }
01194 else if ((bus[tempa].phases & 0x20) == 0x20)
01195 {
01196
01197
01198 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01199 {
01200 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01201 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01202 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yfrom[jindex*3+kindex]).Im());
01203 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01204
01205 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01206 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01207 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(branch[jindexer].Yfrom[jindex*3+kindex]).Im();
01208 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01209 }
01210
01211 if (((branch[jindexer].Yto[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01212 {
01213 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01214 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01215 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Im());
01216 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01217
01218 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01219 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01220 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yto[jindex*3+kindex]).Im();
01221 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01222 }
01223
01224 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01225 {
01226 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01227 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01228 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01229 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01230
01231 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01232 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01233 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01234 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01235 }
01236
01237 if (((branch[jindexer].Yto[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1 && bus[tempb].type != 1))
01238 {
01239 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01240 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01241 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Re());
01242 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01243
01244 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01245 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01246 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Re());
01247 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01248 }
01249 }
01250 else if ((bus[tempb].phases & 0x20) == 0x20)
01251 {
01252
01253
01254
01255 if (branch[jindexer].lnk_type == 1)
01256 {
01257
01258
01259 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01260 {
01261 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01262 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01263 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Im());
01264 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01265
01266 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01267 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01268 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yfrom[jindex*3+kindex]).Im();
01269 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01270 }
01271
01272 if (((branch[jindexer].Yto[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01273 {
01274 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01275 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01276 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yto[jindex*3+kindex]).Im());
01277 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01278
01279 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01280 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01281 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(branch[jindexer].Yto[jindex*3+kindex]).Im();
01282 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01283 }
01284
01285 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01286 {
01287 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01288 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01289 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01290 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01291
01292 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01293 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01294 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01295 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01296 }
01297
01298 if (((branch[jindexer].Yto[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1 && bus[tempb].type != 1))
01299 {
01300 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01301 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01302 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yto[jindex*3+kindex]).Re());
01303 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01304
01305 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01306 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01307 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yto[jindex*3+kindex]).Re());
01308 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01309 }
01310 }
01311 else
01312 {
01313 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01314 {
01315 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01316 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01317 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Im());
01318 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01319
01320 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01321 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01322 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yfrom[jindex*3+kindex]).Im();
01323 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01324 }
01325
01326 if (((branch[jindexer].Yto[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01327 {
01328 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01329 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01330 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yto[jindex*3+kindex]).Im());
01331 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01332
01333 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01334 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01335 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(branch[jindexer].Yto[jindex*3+kindex]).Im();
01336 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01337 }
01338
01339 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01340 {
01341 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01342 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01343 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01344 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01345
01346 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01347 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01348 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01349 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01350 }
01351
01352 if (((branch[jindexer].Yto[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1 && bus[tempb].type != 1))
01353 {
01354 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01355 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01356 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yto[jindex*3+kindex]).Re());
01357 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01358
01359 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01360 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01361 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = ((branch[jindexer].Yto[jindex*3+kindex]).Re());
01362 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01363 }
01364 }
01365 }
01366 else
01367 {
01368
01369
01370 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01371 {
01372 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01373 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01374 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Im());
01375 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01376
01377 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01378 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01379 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yfrom[jindex*3+kindex]).Im();
01380 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01381 }
01382
01383 if (((branch[jindexer].Yto[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01384 {
01385 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01386 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01387 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Im());
01388 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01389
01390 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01391 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01392 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yto[jindex*3+kindex]).Im();
01393 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01394 }
01395
01396 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01397 {
01398 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex + 2;
01399 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01400 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01401 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01402
01403 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + jindex;
01404 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01405 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01406 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01407 }
01408
01409 if (((branch[jindexer].Yto[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1 && bus[tempb].type != 1))
01410 {
01411 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex + 2;
01412 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex;
01413 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Re());
01414 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01415
01416 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + jindex;
01417 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + kindex + 2;
01418 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[jindex*3+kindex]).Re());
01419 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01420 }
01421 }
01422 }
01423 }
01424 }
01425 else if ((bus[tempa].phases & 0x80) == 0x80)
01426 {
01427 GL_THROW("NR does not support triplex to 3-phase connections.");
01428
01429
01430
01431
01432
01433
01434 }
01435 else
01436 {
01437
01438 phase_workc = (branch[jindexer].phases & 0x07);
01439
01440
01441 temp_index = -1;
01442 temp_size = -1;
01443
01444
01445 switch(bus[tempa].phases & 0x07)
01446 {
01447 case 0x01:
01448 {
01449 temp_size = 1;
01450
01451 if (phase_workc==0x01)
01452 {
01453
01454 temp_index = 0;
01455 }
01456 else if (phase_workc==0x02)
01457 {
01458 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01459
01460
01461
01462
01463
01464 }
01465 else
01466 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01467
01468 break;
01469 }
01470 case 0x02:
01471 {
01472 temp_size = 1;
01473
01474 if (phase_workc==0x01)
01475 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01476 else if (phase_workc==0x02)
01477 {
01478
01479 temp_index = 0;
01480 }
01481 else
01482 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01483
01484 break;
01485 }
01486 case 0x03:
01487 {
01488 temp_size = 2;
01489
01490 if (phase_workc==0x01)
01491 {
01492
01493 temp_index = 1;
01494 }
01495 else if (phase_workc==0x02)
01496 {
01497
01498 temp_index = 0;
01499 }
01500 else
01501 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01502
01503 break;
01504 }
01505 case 0x04:
01506 {
01507 temp_size = 1;
01508
01509 if (phase_workc==0x01)
01510 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01511 else if (phase_workc==0x02)
01512 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01513 else
01514 {
01515
01516 temp_index = 0;
01517 }
01518
01519 break;
01520 }
01521 case 0x05:
01522 {
01523 temp_size = 2;
01524
01525 if (phase_workc==0x01)
01526 {
01527
01528 temp_index = 1;
01529 }
01530 else if (phase_workc==0x02)
01531 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01532 else
01533 {
01534
01535 temp_index = 0;
01536 }
01537
01538 break;
01539 }
01540 case 0x06:
01541 {
01542 temp_size = 2;
01543
01544 if (phase_workc==0x01)
01545 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01546 else if (phase_workc==0x02)
01547 {
01548
01549 temp_index = 1;
01550 }
01551 else
01552 {
01553
01554 temp_index = 0;
01555 }
01556
01557 break;
01558 }
01559 case 0x07:
01560 {
01561 temp_size = 3;
01562
01563 if (phase_workc==0x01)
01564 {
01565
01566 temp_index = 2;
01567 }
01568 else if (phase_workc==0x02)
01569 {
01570
01571 temp_index = 1;
01572 }
01573 else
01574 {
01575
01576 temp_index = 0;
01577 }
01578
01579 break;
01580 }
01581 default:
01582 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01583 break;
01584 }
01585 if ((temp_index==-1) || (temp_size==-1))
01586 GL_THROW("NR: A center-tapped transformer has an invalid phase matching");
01587
01588
01589 if (phase_workc==0x01)
01590 {
01591 jindex=2;
01592 }
01593 else if (phase_workc==0x02)
01594 {
01595 jindex=1;
01596 }
01597 else
01598 {
01599 jindex=0;
01600 }
01601
01602
01603
01604 for (kindex=0; kindex<2; kindex++)
01605 {
01606
01607 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01608 {
01609 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index;
01610 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01611 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Im());
01612 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01613
01614 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + temp_size;
01615 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01616 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yfrom[jindex*3+kindex]).Im();
01617 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01618 }
01619
01620 if (((branch[jindexer].Yto[kindex*3+jindex]).Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01621 {
01622 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + kindex;
01623 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index;
01624 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[kindex*3+jindex]).Im());
01625 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01626
01627 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01628 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + temp_size;
01629 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (branch[jindexer].Yto[kindex*3+jindex]).Im();
01630 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01631 }
01632
01633 if (((branch[jindexer].Yfrom[jindex*3+kindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01634 {
01635 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + temp_size;
01636 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex;
01637 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01638 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01639
01640 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index;
01641 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01642 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yfrom[jindex*3+kindex]).Re());
01643 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01644 }
01645
01646 if (((branch[jindexer].Yto[kindex*3+jindex]).Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
01647 {
01648 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + kindex + 2;
01649 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index;
01650 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[kindex*3+jindex]).Re());
01651 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01652
01653 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + kindex;
01654 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + temp_size;
01655 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -((branch[jindexer].Yto[kindex*3+jindex]).Re());
01656 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
01657 }
01658 }
01659
01660 }
01661 }
01662 else
01663 {
01664
01665 temp_index = temp_index_b = -1;
01666 temp_size = temp_size_b = temp_size_c = -1;
01667 Full_Mat_A = Full_Mat_B = false;
01668
01669
01670 switch(branch[jindexer].phases & 0x07) {
01671 case 0x00:
01672 {
01673 temp_size_c = -99;
01674 break;
01675 }
01676 case 0x01:
01677 {
01678 Temp_Ad_A[0][0] = branch[jindexer].Yfrom[8];
01679 Temp_Ad_B[0][0] = branch[jindexer].Yto[8];
01680 temp_size_c = 1;
01681 break;
01682 }
01683 case 0x02:
01684 {
01685 Temp_Ad_A[0][0] = branch[jindexer].Yfrom[4];
01686 Temp_Ad_B[0][0] = branch[jindexer].Yto[4];
01687 temp_size_c = 1;
01688 break;
01689 }
01690 case 0x03:
01691 {
01692 Temp_Ad_A[0][0] = branch[jindexer].Yfrom[4];
01693 Temp_Ad_A[0][1] = branch[jindexer].Yfrom[5];
01694 Temp_Ad_A[1][0] = branch[jindexer].Yfrom[7];
01695 Temp_Ad_A[1][1] = branch[jindexer].Yfrom[8];
01696
01697 Temp_Ad_B[0][0] = branch[jindexer].Yto[4];
01698 Temp_Ad_B[0][1] = branch[jindexer].Yto[5];
01699 Temp_Ad_B[1][0] = branch[jindexer].Yto[7];
01700 Temp_Ad_B[1][1] = branch[jindexer].Yto[8];
01701
01702 temp_size_c = 2;
01703 break;
01704 }
01705 case 0x04:
01706 {
01707 Temp_Ad_A[0][0] = branch[jindexer].Yfrom[0];
01708 Temp_Ad_B[0][0] = branch[jindexer].Yto[0];
01709 temp_size_c = 1;
01710 break;
01711 }
01712 case 0x05:
01713 {
01714 Temp_Ad_A[0][0] = branch[jindexer].Yfrom[0];
01715 Temp_Ad_A[0][1] = branch[jindexer].Yfrom[2];
01716 Temp_Ad_A[1][0] = branch[jindexer].Yfrom[6];
01717 Temp_Ad_A[1][1] = branch[jindexer].Yfrom[8];
01718
01719 Temp_Ad_B[0][0] = branch[jindexer].Yto[0];
01720 Temp_Ad_B[0][1] = branch[jindexer].Yto[2];
01721 Temp_Ad_B[1][0] = branch[jindexer].Yto[6];
01722 Temp_Ad_B[1][1] = branch[jindexer].Yto[8];
01723
01724 temp_size_c = 2;
01725 break;
01726 }
01727 case 0x06:
01728 {
01729 Temp_Ad_A[0][0] = branch[jindexer].Yfrom[0];
01730 Temp_Ad_A[0][1] = branch[jindexer].Yfrom[1];
01731 Temp_Ad_A[1][0] = branch[jindexer].Yfrom[3];
01732 Temp_Ad_A[1][1] = branch[jindexer].Yfrom[4];
01733
01734 Temp_Ad_B[0][0] = branch[jindexer].Yto[0];
01735 Temp_Ad_B[0][1] = branch[jindexer].Yto[1];
01736 Temp_Ad_B[1][0] = branch[jindexer].Yto[3];
01737 Temp_Ad_B[1][1] = branch[jindexer].Yto[4];
01738
01739 temp_size_c = 2;
01740 break;
01741 }
01742 default:
01743 {
01744 break;
01745 }
01746 }
01747
01748 if (temp_size_c==-99)
01749 {
01750 continue;
01751 }
01752
01753 if (temp_size_c==-1)
01754 {
01755 GL_THROW("NR: A line's phase was flagged as not full three-phase, but wasn't: (%s) %u %u",
01756 branch[jindexer].name, branch[jindexer].phases, branch[jindexer].origphases);
01757
01758
01759
01760
01761
01762 }
01763
01764
01765 switch(bus[tempa].phases & 0x07) {
01766 case 0x01:
01767 {
01768 if ((branch[jindexer].phases & 0x07) == 0x01)
01769 {
01770 temp_size = 1;
01771 temp_index = 0;
01772 }
01773 else
01774 {
01775 GL_THROW("NR: One of the lines has invalid phase parameters");
01776
01777
01778
01779
01780
01781 }
01782 break;
01783 }
01784 case 0x02:
01785 {
01786 if ((branch[jindexer].phases & 0x07) == 0x02)
01787 {
01788 temp_size = 1;
01789 temp_index = 0;
01790 }
01791 else
01792 {
01793 GL_THROW("NR: One of the lines has invalid phase parameters");
01794 }
01795 break;
01796 }
01797 case 0x03:
01798 {
01799 temp_size = 2;
01800 if ((branch[jindexer].phases & 0x07) == 0x01)
01801 {
01802 temp_index = 1;
01803 }
01804 else if ((branch[jindexer].phases & 0x07) == 0x02)
01805 {
01806 temp_index = 0;
01807 }
01808 else if ((branch[jindexer].phases & 0x07) == 0x03)
01809 {
01810 temp_index = 0;
01811 }
01812 else
01813 {
01814 GL_THROW("NR: One of the lines has invalid phase parameters");
01815 }
01816 break;
01817 }
01818 case 0x04:
01819 {
01820 if ((branch[jindexer].phases & 0x07) == 0x04)
01821 {
01822 temp_size = 1;
01823 temp_index = 0;
01824 }
01825 else
01826 {
01827 GL_THROW("NR: One of the lines has invalid phase parameters");
01828 }
01829 break;
01830 }
01831 case 0x05:
01832 {
01833 temp_size = 2;
01834 if ((branch[jindexer].phases & 0x07) == 0x01)
01835 {
01836 temp_index = 1;
01837 }
01838 else if ((branch[jindexer].phases & 0x07) == 0x04)
01839 {
01840 temp_index = 0;
01841 }
01842 else if ((branch[jindexer].phases & 0x07) == 0x05)
01843 {
01844 temp_index = 0;
01845 }
01846 else
01847 {
01848 GL_THROW("NR: One of the lines has invalid phase parameters");
01849 }
01850 break;
01851 }
01852 case 0x06:
01853 {
01854 temp_size = 2;
01855 if ((branch[jindexer].phases & 0x07) == 0x02)
01856 {
01857 temp_index = 1;
01858 }
01859 else if ((branch[jindexer].phases & 0x07) == 0x04)
01860 {
01861 temp_index = 0;
01862 }
01863 else if ((branch[jindexer].phases & 0x07) == 0x06)
01864 {
01865 temp_index = 0;
01866 }
01867 else
01868 {
01869 GL_THROW("NR: One of the lines has invalid phase parameters");
01870 }
01871 break;
01872 }
01873 case 0x07:
01874 {
01875 temp_size = 3;
01876 if ((branch[jindexer].phases & 0x07) == 0x01)
01877 {
01878 temp_index = 2;
01879 }
01880 else if ((branch[jindexer].phases & 0x07) == 0x02)
01881 {
01882 temp_index = 1;
01883 }
01884 else if ((branch[jindexer].phases & 0x07) == 0x03)
01885 {
01886 temp_index = 1;
01887 }
01888 else if ((branch[jindexer].phases & 0x07) == 0x04)
01889 {
01890 temp_index = 0;
01891 }
01892 else if ((branch[jindexer].phases & 0x07) == 0x05)
01893 {
01894 temp_index = 0;
01895 Full_Mat_A = true;
01896 }
01897 else if ((branch[jindexer].phases & 0x07) == 0x06)
01898 {
01899 temp_index = 0;
01900 }
01901 else
01902 {
01903 GL_THROW("NR: One of the lines has invalid phase parameters");
01904 }
01905 break;
01906 }
01907 default:
01908 {
01909 break;
01910 }
01911 }
01912
01913
01914 switch(bus[tempb].phases & 0x07) {
01915 case 0x01:
01916 {
01917 if ((branch[jindexer].phases & 0x07) == 0x01)
01918 {
01919 temp_size_b = 1;
01920 temp_index_b = 0;
01921 }
01922 else
01923 {
01924 GL_THROW("NR: One of the lines has invalid phase parameters");
01925
01926
01927
01928
01929
01930 }
01931 break;
01932 }
01933 case 0x02:
01934 {
01935 if ((branch[jindexer].phases & 0x07) == 0x02)
01936 {
01937 temp_size_b = 1;
01938 temp_index_b = 0;
01939 }
01940 else
01941 {
01942 GL_THROW("NR: One of the lines has invalid phase parameters");
01943 }
01944 break;
01945 }
01946 case 0x03:
01947 {
01948 temp_size_b = 2;
01949 if ((branch[jindexer].phases & 0x07) == 0x01)
01950 {
01951 temp_index_b = 1;
01952 }
01953 else if ((branch[jindexer].phases & 0x07) == 0x02)
01954 {
01955 temp_index_b = 0;
01956 }
01957 else if ((branch[jindexer].phases & 0x07) == 0x03)
01958 {
01959 temp_index_b = 0;
01960 }
01961 else
01962 {
01963 GL_THROW("NR: One of the lines has invalid phase parameters");
01964 }
01965 break;
01966 }
01967 case 0x04:
01968 {
01969 if ((branch[jindexer].phases & 0x07) == 0x04)
01970 {
01971 temp_size_b = 1;
01972 temp_index_b = 0;
01973 }
01974 else
01975 {
01976 GL_THROW("NR: One of the lines has invalid phase parameters");
01977 }
01978 break;
01979 }
01980 case 0x05:
01981 {
01982 temp_size_b = 2;
01983 if ((branch[jindexer].phases & 0x07) == 0x01)
01984 {
01985 temp_index_b = 1;
01986 }
01987 else if ((branch[jindexer].phases & 0x07) == 0x04)
01988 {
01989 temp_index_b = 0;
01990 }
01991 else if ((branch[jindexer].phases & 0x07) == 0x05)
01992 {
01993 temp_index_b = 0;
01994 }
01995 else
01996 {
01997 GL_THROW("NR: One of the lines has invalid phase parameters");
01998 }
01999 break;
02000 }
02001 case 0x06:
02002 {
02003 temp_size_b = 2;
02004 if ((branch[jindexer].phases & 0x07) == 0x02)
02005 {
02006 temp_index_b = 1;
02007 }
02008 else if ((branch[jindexer].phases & 0x07) == 0x04)
02009 {
02010 temp_index_b = 0;
02011 }
02012 else if ((branch[jindexer].phases & 0x07) == 0x06)
02013 {
02014 temp_index_b = 0;
02015 }
02016 else
02017 {
02018 GL_THROW("NR: One of the lines has invalid phase parameters");
02019 }
02020 break;
02021 }
02022 case 0x07:
02023 {
02024 temp_size_b = 3;
02025 if ((branch[jindexer].phases & 0x07) == 0x01)
02026 {
02027 temp_index_b = 2;
02028 }
02029 else if ((branch[jindexer].phases & 0x07) == 0x02)
02030 {
02031 temp_index_b = 1;
02032 }
02033 else if ((branch[jindexer].phases & 0x07) == 0x03)
02034 {
02035 temp_index_b = 1;
02036 }
02037 else if ((branch[jindexer].phases & 0x07) == 0x04)
02038 {
02039 temp_index_b = 0;
02040 }
02041 else if ((branch[jindexer].phases & 0x07) == 0x05)
02042 {
02043 temp_index_b = 0;
02044 Full_Mat_B = true;
02045 }
02046 else if ((branch[jindexer].phases & 0x07) == 0x06)
02047 {
02048 temp_index_b = 0;
02049 }
02050 else
02051 {
02052 GL_THROW("NR: One of the lines has invalid phase parameters");
02053 }
02054 break;
02055 }
02056 default:
02057 {
02058 break;
02059 }
02060 }
02061
02062
02063 if ((temp_index==-1) || (temp_index_b==-1) || (temp_size==-1) || (temp_size_b==-1) || (temp_size_c==-1))
02064 {
02065 GL_THROW("NR: Failure to construct single/double phase line indices");
02066
02067
02068
02069
02070 }
02071
02072 if (Full_Mat_A)
02073 {
02074 for (jindex=0; jindex<temp_size_c; jindex++)
02075 {
02076 for (kindex=0; kindex<temp_size_c; kindex++)
02077 {
02078
02079 if ((Temp_Ad_A[jindex][kindex].Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02080 {
02081 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex*2;
02082 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex;
02083 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Im());
02084 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02085
02086 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex*2 + temp_size;
02087 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex + temp_size_b;
02088 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (Temp_Ad_A[jindex][kindex].Im());
02089 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02090 }
02091
02092 if ((Temp_Ad_B[jindex][kindex].Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02093 {
02094 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex;
02095 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex*2;
02096 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Im());
02097 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02098
02099 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex + temp_size_b;
02100 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex*2 + temp_size;
02101 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = Temp_Ad_B[jindex][kindex].Im();
02102 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02103 }
02104
02105 if ((Temp_Ad_A[jindex][kindex].Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02106 {
02107 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex*2 + temp_size;
02108 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex;
02109 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Re());
02110 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02111
02112 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex*2;
02113 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex + temp_size_b;
02114 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Re());
02115 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02116 }
02117
02118 if ((Temp_Ad_B[jindex][kindex].Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02119 {
02120 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex + temp_size_b;
02121 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex*2;
02122 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Re());
02123 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02124
02125 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex;
02126 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex*2 + temp_size;
02127 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Re());
02128 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02129 }
02130 }
02131 }
02132 }
02133
02134 if (Full_Mat_B)
02135 {
02136 for (jindex=0; jindex<temp_size_c; jindex++)
02137 {
02138 for (kindex=0; kindex<temp_size_c; kindex++)
02139 {
02140
02141 if ((Temp_Ad_A[jindex][kindex].Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02142 {
02143 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex;
02144 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex*2;
02145 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Im());
02146 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02147
02148 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex + temp_size;
02149 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex*2 + temp_size_b;
02150 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (Temp_Ad_A[jindex][kindex].Im());
02151 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02152 }
02153
02154 if ((Temp_Ad_B[jindex][kindex].Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02155 {
02156 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex*2;
02157 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex;
02158 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Im());
02159 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02160
02161 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex*2 + temp_size_b;
02162 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex + temp_size;
02163 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = Temp_Ad_B[jindex][kindex].Im();
02164 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02165 }
02166
02167 if ((Temp_Ad_A[jindex][kindex].Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02168 {
02169 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex + temp_size;
02170 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex*2;
02171 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Re());
02172 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02173
02174 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex;
02175 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex*2 + temp_size_b;
02176 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Re());
02177 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02178 }
02179
02180 if ((Temp_Ad_B[jindex][kindex].Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02181 {
02182 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex*2 + temp_size_b;
02183 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex;
02184 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Re());
02185 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02186
02187 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex*2;
02188 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex + temp_size;
02189 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Re());
02190 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02191 }
02192 }
02193 }
02194 }
02195
02196 if ((!Full_Mat_A) && (!Full_Mat_B))
02197 {
02198 for (jindex=0; jindex<temp_size_c; jindex++)
02199 {
02200 for (kindex=0; kindex<temp_size_c; kindex++)
02201 {
02202
02203 if ((Temp_Ad_A[jindex][kindex].Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02204 {
02205 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex;
02206 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex;
02207 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Im());
02208 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02209
02210 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex + temp_size;
02211 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex + temp_size_b;
02212 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (Temp_Ad_A[jindex][kindex].Im());
02213 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02214 }
02215
02216 if ((Temp_Ad_B[jindex][kindex].Im() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02217 {
02218 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex;
02219 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex;
02220 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Im());
02221 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02222
02223 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex + temp_size_b;
02224 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex + temp_size;
02225 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = Temp_Ad_B[jindex][kindex].Im();
02226 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02227 }
02228
02229 if ((Temp_Ad_A[jindex][kindex].Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02230 {
02231 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex + temp_size;
02232 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex;
02233 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Re());
02234 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02235
02236 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempa].Matrix_Loc + temp_index + jindex;
02237 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + kindex + temp_size_b;
02238 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_A[jindex][kindex].Re());
02239 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02240 }
02241
02242 if ((Temp_Ad_B[jindex][kindex].Re() != 0) && (bus[tempa].type != 1) && (bus[tempb].type != 1))
02243 {
02244 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex + temp_size_b;
02245 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex;
02246 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Re());
02247 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02248
02249 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*bus[tempb].Matrix_Loc + temp_index_b + jindex;
02250 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*bus[tempa].Matrix_Loc + temp_index + kindex + temp_size;
02251 powerflow_values->island_matrix_values[island_index_val].Y_offdiag_PQ[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(Temp_Ad_B[jindex][kindex].Re());
02252 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02253 }
02254 }
02255 }
02256 }
02257 }
02258 }
02259
02260
02261 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
02262 {
02263
02264 powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed = 0;
02265 }
02266
02267 for (jindexer=0; jindexer<bus_count;jindexer++)
02268 {
02269 if (bus[jindexer].island_number != -1)
02270 {
02271
02272 island_loop_index = bus[jindexer].island_number;
02273
02274
02275 for (jindex=0; jindex<3; jindex++)
02276 {
02277 for (kindex=0; kindex<3; kindex++)
02278 {
02279 if ((powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Re() != 0 && bus[jindexer].type != 1 && jindex!=kindex)
02280 powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed += 1;
02281 if ((powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Im() != 0 && bus[jindexer].type != 1 && jindex!=kindex)
02282 powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed += 1;
02283 }
02284 }
02285 }
02286
02287 }
02288
02289
02290 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
02291 {
02292 if (powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed == NULL)
02293 {
02294 powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed = (Y_NR *)gl_malloc((powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed*2) *sizeof(Y_NR));
02295
02296
02297 if (powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed == NULL)
02298 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
02299
02300
02301 powerflow_values->island_matrix_values[island_loop_index].max_size_diag_fixed = powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed;
02302 }
02303 else if (powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed > powerflow_values->island_matrix_values[island_loop_index].max_size_diag_fixed)
02304 {
02305
02306 gl_free(powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed);
02307
02308
02309 powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed = (Y_NR *)gl_malloc((powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed*2) *sizeof(Y_NR));
02310
02311
02312 if (powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed == NULL)
02313 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
02314
02315
02316 powerflow_values->island_matrix_values[island_loop_index].max_size_diag_fixed = powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed;
02317
02318
02319 powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed = true;
02320 }
02321
02322
02323 powerflow_values->island_matrix_values[island_loop_index].indexer = 0;
02324 }
02325
02326 for (jindexer=0; jindexer<bus_count;jindexer++)
02327 {
02328
02329 island_index_val = bus[jindexer].island_number;
02330
02331
02332 if (island_index_val == -1)
02333 {
02334 continue;
02335 }
02336
02337
02338 for (jindex=0; jindex<powerflow_values->BA_diag[jindexer].size; jindex++)
02339 {
02340 for (kindex=0; kindex<powerflow_values->BA_diag[jindexer].size; kindex++)
02341 {
02342 if ((powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Im() != 0 && bus[jindexer].type != 1 && jindex!=kindex)
02343 {
02344 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*powerflow_values->BA_diag[jindexer].row_ind + jindex;
02345 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*powerflow_values->BA_diag[jindexer].col_ind + kindex;
02346 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Im();
02347 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02348
02349 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*powerflow_values->BA_diag[jindexer].row_ind + jindex +powerflow_values->BA_diag[jindexer].size;
02350 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*powerflow_values->BA_diag[jindexer].col_ind + kindex +powerflow_values->BA_diag[jindexer].size;
02351 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = -(powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Im();
02352 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02353 }
02354
02355 if ((powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Re() != 0 && bus[jindexer].type != 1 && jindex!=kindex)
02356 {
02357 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*powerflow_values->BA_diag[jindexer].row_ind + jindex;
02358 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*powerflow_values->BA_diag[jindexer].col_ind + kindex +powerflow_values->BA_diag[jindexer].size;
02359 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Re();
02360 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02361
02362 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].row_ind = 2*powerflow_values->BA_diag[jindexer].row_ind + jindex +powerflow_values->BA_diag[jindexer].size;
02363 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].col_ind = 2*powerflow_values->BA_diag[jindexer].col_ind + kindex;
02364 powerflow_values->island_matrix_values[island_index_val].Y_diag_fixed[powerflow_values->island_matrix_values[island_index_val].indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][kindex]).Re();
02365 powerflow_values->island_matrix_values[island_index_val].indexer += 1;
02366 }
02367 }
02368 }
02369 }
02370 }
02371
02372
02373 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
02374 {
02375
02376 powerflow_values->island_matrix_values[island_loop_index].SaturationMismatchPresent = false;
02377
02378
02379 powerflow_values->island_matrix_values[island_loop_index].iteration_count = 0;
02380 }
02381
02382
02383 still_iterating_islands = true;
02384
02385
02386 island_loop_index = 0;
02387
02388
02389 while (still_iterating_islands == true)
02390 {
02391
02392 if (matrix_solver_method==MM_SUPERLU)
02393 {
02394
02395 curr_island_superLU_vars = (SUPERLU_NR_vars *)powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars;
02396 }
02397
02398
02399 compute_load_values(bus_count,bus,powerflow_values,false,island_loop_index);
02400
02401
02402
02403 if (powerflow_values->island_matrix_values[island_loop_index].deltaI_NR==NULL)
02404 {
02405 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR = (double *)gl_malloc((2*powerflow_values->island_matrix_values[island_loop_index].total_variables) *sizeof(double));
02406
02407
02408 if (powerflow_values->island_matrix_values[island_loop_index].deltaI_NR == NULL)
02409 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
02410
02411
02412 powerflow_values->island_matrix_values[island_loop_index].max_total_variables = powerflow_values->island_matrix_values[island_loop_index].total_variables;
02413 }
02414 else if (powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed)
02415 {
02416
02417 gl_free(powerflow_values->island_matrix_values[island_loop_index].deltaI_NR);
02418
02419
02420 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR = (double *)gl_malloc((2*powerflow_values->island_matrix_values[island_loop_index].total_variables) *sizeof(double));
02421
02422
02423 if (powerflow_values->island_matrix_values[island_loop_index].deltaI_NR == NULL)
02424 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
02425
02426
02427 powerflow_values->island_matrix_values[island_loop_index].max_total_variables = powerflow_values->island_matrix_values[island_loop_index].total_variables;
02428 }
02429
02430
02431 if (powerflow_type == PF_DYNINIT)
02432 {
02433 powerflow_values->island_matrix_values[island_loop_index].swing_converged = true;
02434 }
02435
02436
02437 for (indexer=0; indexer<bus_count; indexer++)
02438 {
02439
02440 if (bus[indexer].island_number == island_loop_index)
02441 {
02442
02443 if ((powerflow_type == PF_DYNINIT) && (powerflow_values->island_matrix_values[island_loop_index].swing_is_a_swing==true))
02444 {
02445
02446
02447 if ((*bus[indexer].dynamics_enabled==true) && (bus[indexer].full_Y != NULL) && (bus[indexer].DynCurrent != NULL))
02448 {
02449
02450 temp_complex_1 = (~bus[indexer].V[0]) + (~bus[indexer].V[1])*avalsq + (~bus[indexer].V[2])*aval;
02451
02452
02453
02454 temp_complex_2 = ~bus[indexer].V[0];
02455
02456
02457 temp_complex_0 = temp_complex_2*(bus[indexer].full_Y[0]*bus[indexer].V[0] + bus[indexer].full_Y[1]*bus[indexer].V[1] + bus[indexer].full_Y[2]*bus[indexer].V[2]);
02458
02459
02460 temp_complex_5 = bus[indexer].full_Y[0]*bus[indexer].V[0] + bus[indexer].full_Y[1]*bus[indexer].V[1] + bus[indexer].full_Y[2]*bus[indexer].V[2];
02461 temp_complex_4 = ~temp_complex_5;
02462 temp_complex_3 = bus[indexer].V[0]*temp_complex_4;
02463
02464
02465 temp_complex_2 = ~bus[indexer].V[1];
02466
02467
02468 temp_complex_0 += temp_complex_2*(bus[indexer].full_Y[3]*bus[indexer].V[0] + bus[indexer].full_Y[4]*bus[indexer].V[1] + bus[indexer].full_Y[5]*bus[indexer].V[2]);
02469
02470
02471 temp_complex_5 = bus[indexer].full_Y[3]*bus[indexer].V[0] + bus[indexer].full_Y[4]*bus[indexer].V[1] + bus[indexer].full_Y[5]*bus[indexer].V[2];
02472 temp_complex_4 = ~temp_complex_5;
02473 temp_complex_3 += bus[indexer].V[1]*temp_complex_4;
02474
02475
02476 temp_complex_2 = ~bus[indexer].V[2];
02477
02478
02479 temp_complex_0 += temp_complex_2*(bus[indexer].full_Y[6]*bus[indexer].V[0] + bus[indexer].full_Y[7]*bus[indexer].V[1] + bus[indexer].full_Y[8]*bus[indexer].V[2]);
02480
02481
02482 temp_complex_5 = bus[indexer].full_Y[6]*bus[indexer].V[0] + bus[indexer].full_Y[7]*bus[indexer].V[1] + bus[indexer].full_Y[8]*bus[indexer].V[2];
02483 temp_complex_4 = ~temp_complex_5;
02484 temp_complex_3 += bus[indexer].V[2]*temp_complex_4;
02485
02486
02487
02488
02489 if ((bus[indexer].type == 2) && (powerflow_values->island_matrix_values[island_loop_index].iteration_count>0))
02490 {
02491 *bus[indexer].PGenTotal = complex(0.0,0.0);
02492 }
02493 }
02494 else
02495 {
02496 temp_complex_0 = 0.0;
02497 temp_complex_1 = 0.0;
02498 }
02499 }
02500
02501 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
02502 {
02503 tempIcalcReal = tempIcalcImag = 0;
02504
02505 if ((bus[indexer].phases & 0x80) == 0x80)
02506 {
02507
02508
02509 if ((bus[indexer].phases & 0x20) == 0x20)
02510 {
02511
02512 tempPbus = bus[indexer].PL[jindex];
02513 tempQbus = bus[indexer].QL[jindex];
02514 }
02515 else
02516 {
02517
02518 tempPbus = -bus[indexer].PL[jindex];
02519 tempQbus = -bus[indexer].QL[jindex];
02520 }
02521
02522
02523
02524 tempIcalcReal += (powerflow_values->BA_diag[indexer].Y[jindex][0]).Re() * (bus[indexer].V[0]).Re() - (powerflow_values->BA_diag[indexer].Y[jindex][0]).Im() * (bus[indexer].V[0]).Im();
02525 tempIcalcImag += (powerflow_values->BA_diag[indexer].Y[jindex][0]).Re() * (bus[indexer].V[0]).Im() + (powerflow_values->BA_diag[indexer].Y[jindex][0]).Im() * (bus[indexer].V[0]).Re();
02526
02527
02528 tempIcalcReal += (powerflow_values->BA_diag[indexer].Y[jindex][1]).Re() * (bus[indexer].V[1]).Re() - (powerflow_values->BA_diag[indexer].Y[jindex][1]).Im() * (bus[indexer].V[1]).Im();
02529 tempIcalcImag += (powerflow_values->BA_diag[indexer].Y[jindex][1]).Re() * (bus[indexer].V[1]).Im() + (powerflow_values->BA_diag[indexer].Y[jindex][1]).Im() * (bus[indexer].V[1]).Re();
02530
02531
02532 for (kindexer=0; kindexer<(bus[indexer].Link_Table_Size); kindexer++)
02533 {
02534
02535 jindexer=bus[indexer].Link_Table[kindexer];
02536
02537 if (branch[jindexer].from == indexer)
02538 {
02539 if ((bus[indexer].phases & 0x20) == 0x20)
02540 {
02541 work_vals_char_0 = jindex*3;
02542
02543
02544
02545
02546 tempIcalcReal += ((branch[jindexer].Yfrom[work_vals_char_0])).Re() * (bus[branch[jindexer].to].V[0]).Re() - ((branch[jindexer].Yfrom[work_vals_char_0])).Im() * (bus[branch[jindexer].to].V[0]).Im();
02547 tempIcalcImag += ((branch[jindexer].Yfrom[work_vals_char_0])).Re() * (bus[branch[jindexer].to].V[0]).Im() + ((branch[jindexer].Yfrom[work_vals_char_0])).Im() * (bus[branch[jindexer].to].V[0]).Re();
02548
02549
02550 tempIcalcReal += ((branch[jindexer].Yfrom[jindex*3+1])).Re() * (bus[branch[jindexer].to].V[1]).Re() - ((branch[jindexer].Yfrom[jindex*3+1])).Im() * (bus[branch[jindexer].to].V[1]).Im();
02551 tempIcalcImag += ((branch[jindexer].Yfrom[jindex*3+1])).Re() * (bus[branch[jindexer].to].V[1]).Im() + ((branch[jindexer].Yfrom[jindex*3+1])).Im() * (bus[branch[jindexer].to].V[1]).Re();
02552
02553 }
02554 else
02555 {
02556 work_vals_char_0 = jindex*3;
02557
02558
02559
02560 tempIcalcReal += (-(branch[jindexer].Yfrom[work_vals_char_0])).Re() * (bus[branch[jindexer].to].V[0]).Re() - (-(branch[jindexer].Yfrom[work_vals_char_0])).Im() * (bus[branch[jindexer].to].V[0]).Im();
02561 tempIcalcImag += (-(branch[jindexer].Yfrom[work_vals_char_0])).Re() * (bus[branch[jindexer].to].V[0]).Im() + (-(branch[jindexer].Yfrom[work_vals_char_0])).Im() * (bus[branch[jindexer].to].V[0]).Re();
02562
02563
02564 tempIcalcReal += (-(branch[jindexer].Yfrom[jindex*3+1])).Re() * (bus[branch[jindexer].to].V[1]).Re() - (-(branch[jindexer].Yfrom[jindex*3+1])).Im() * (bus[branch[jindexer].to].V[1]).Im();
02565 tempIcalcImag += (-(branch[jindexer].Yfrom[jindex*3+1])).Re() * (bus[branch[jindexer].to].V[1]).Im() + (-(branch[jindexer].Yfrom[jindex*3+1])).Im() * (bus[branch[jindexer].to].V[1]).Re();
02566
02567 }
02568 }
02569 else if (branch[jindexer].to == indexer)
02570 {
02571 if (branch[jindexer].v_ratio != 1.0)
02572 {
02573
02574 if ((branch[jindexer].phases & 0x01) == 0x01)
02575 {
02576 temp_index=2;
02577 }
02578 else if ((branch[jindexer].phases & 0x02) == 0x02)
02579 {
02580 temp_index=1;
02581 }
02582 else if ((branch[jindexer].phases & 0x04) == 0x04)
02583 {
02584 temp_index=0;
02585 }
02586 else
02587 {
02588 GL_THROW("NR: A split-phase transformer appears to have an invalid phase");
02589 }
02590
02591 work_vals_char_0 = jindex*3+temp_index;
02592
02593
02594 tempIcalcReal += (-(branch[jindexer].Yto[work_vals_char_0])).Re() * (bus[branch[jindexer].from].V[temp_index]).Re() - (-(branch[jindexer].Yto[work_vals_char_0])).Im() * (bus[branch[jindexer].from].V[temp_index]).Im();
02595 tempIcalcImag += (-(branch[jindexer].Yto[work_vals_char_0])).Re() * (bus[branch[jindexer].from].V[temp_index]).Im() + (-(branch[jindexer].Yto[work_vals_char_0])).Im() * (bus[branch[jindexer].from].V[temp_index]).Re();
02596
02597 }
02598 else
02599 {
02600 if ((bus[indexer].phases & 0x20) == 0x20)
02601 {
02602 work_vals_char_0 = jindex*3;
02603
02604
02605
02606 tempIcalcReal += ((branch[jindexer].Yto[work_vals_char_0])).Re() * (bus[branch[jindexer].from].V[0]).Re() - ((branch[jindexer].Yto[work_vals_char_0])).Im() * (bus[branch[jindexer].from].V[0]).Im();
02607 tempIcalcImag += ((branch[jindexer].Yto[work_vals_char_0])).Re() * (bus[branch[jindexer].from].V[0]).Im() + ((branch[jindexer].Yto[work_vals_char_0])).Im() * (bus[branch[jindexer].from].V[0]).Re();
02608
02609
02610 tempIcalcReal += ((branch[jindexer].Yto[work_vals_char_0+1])).Re() * (bus[branch[jindexer].from].V[1]).Re() - ((branch[jindexer].Yto[work_vals_char_0+1])).Im() * (bus[branch[jindexer].from].V[1]).Im();
02611 tempIcalcImag += ((branch[jindexer].Yto[work_vals_char_0+1])).Re() * (bus[branch[jindexer].from].V[1]).Im() + ((branch[jindexer].Yto[work_vals_char_0+1])).Im() * (bus[branch[jindexer].from].V[1]).Re();
02612 }
02613 else
02614 {
02615 work_vals_char_0 = jindex*3;
02616
02617
02618 tempIcalcReal += (-(branch[jindexer].Yto[work_vals_char_0])).Re() * (bus[branch[jindexer].from].V[0]).Re() - (-(branch[jindexer].Yto[work_vals_char_0])).Im() * (bus[branch[jindexer].from].V[0]).Im();
02619 tempIcalcImag += (-(branch[jindexer].Yto[work_vals_char_0])).Re() * (bus[branch[jindexer].from].V[0]).Im() + (-(branch[jindexer].Yto[work_vals_char_0])).Im() * (bus[branch[jindexer].from].V[0]).Re();
02620
02621
02622 tempIcalcReal += (-(branch[jindexer].Yto[work_vals_char_0+1])).Re() * (bus[branch[jindexer].from].V[1]).Re() - (-(branch[jindexer].Yto[work_vals_char_0+1])).Im() * (bus[branch[jindexer].from].V[1]).Im();
02623 tempIcalcImag += (-(branch[jindexer].Yto[work_vals_char_0+1])).Re() * (bus[branch[jindexer].from].V[1]).Im() + (-(branch[jindexer].Yto[work_vals_char_0+1])).Im() * (bus[branch[jindexer].from].V[1]).Re();
02624 }
02625 }
02626 }
02627 else
02628 ;
02629
02630 }
02631
02632
02633 if ((bus[indexer].type > 1) && (bus[indexer].swing_functions_enabled == true))
02634 {
02635 if ((powerflow_type == PF_DYNINIT) && (bus[indexer].DynCurrent != NULL))
02636 {
02637
02638 work_vals_double_3 = tempPbus - tempIcalcReal;
02639 work_vals_double_4 = tempQbus - tempIcalcImag;
02640
02641
02642 work_vals_double_0 = sqrt((work_vals_double_3*work_vals_double_3+work_vals_double_4*work_vals_double_4));
02643
02644 if (work_vals_double_0 > bus[indexer].max_volt_error)
02645 {
02646 powerflow_values->island_matrix_values[island_loop_index].swing_converged=false;
02647 }
02648 }
02649
02650
02651
02652 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + jindex] = 0.0;
02653 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = 0.0;
02654 }
02655 else
02656 {
02657
02658 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+ powerflow_values->BA_diag[indexer].size + jindex] = tempPbus - tempIcalcReal;
02659 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = tempQbus - tempIcalcImag;
02660 }
02661 }
02662 else
02663 {
02664 tempPbus = - bus[indexer].PL[jindex];
02665 tempQbus = - bus[indexer].QL[jindex];
02666
02667 for (kindex=0; kindex<powerflow_values->BA_diag[indexer].size; kindex++)
02668 {
02669
02670 temp_index = -1;
02671 switch(bus[indexer].phases & 0x07) {
02672 case 0x01:
02673 {
02674 temp_index=2;
02675 break;
02676 }
02677 case 0x02:
02678 {
02679 temp_index=1;
02680 break;
02681 }
02682 case 0x03:
02683 {
02684 if (kindex==0)
02685 temp_index=1;
02686 else
02687 temp_index=2;
02688 break;
02689 }
02690 case 0x04:
02691 {
02692 temp_index=0;
02693 break;
02694 }
02695 case 0x05:
02696 {
02697 if (kindex==0)
02698 temp_index=0;
02699 else
02700 temp_index=2;
02701 break;
02702 }
02703 case 0x06:
02704 {
02705 if (kindex==0)
02706 temp_index=0;
02707 else
02708 temp_index=1;
02709 break;
02710 }
02711 case 0x07:
02712 {
02713 temp_index = kindex;
02714 break;
02715 }
02716 }
02717
02718 if (temp_index==-1)
02719 {
02720 GL_THROW("NR: A voltage index failed to be found.");
02721
02722
02723
02724
02725 }
02726
02727
02728 tempIcalcReal += (powerflow_values->BA_diag[indexer].Y[jindex][kindex]).Re() * (bus[indexer].V[temp_index]).Re() - (powerflow_values->BA_diag[indexer].Y[jindex][kindex]).Im() * (bus[indexer].V[temp_index]).Im();
02729 tempIcalcImag += (powerflow_values->BA_diag[indexer].Y[jindex][kindex]).Re() * (bus[indexer].V[temp_index]).Im() + (powerflow_values->BA_diag[indexer].Y[jindex][kindex]).Im() * (bus[indexer].V[temp_index]).Re();
02730
02731
02732 if ((bus[indexer].full_Y_load != NULL) && (jindex==kindex))
02733 {
02734 tempIcalcReal += (bus[indexer].full_Y_load[temp_index]).Re() * (bus[indexer].V[temp_index]).Re() - (bus[indexer].full_Y_load[temp_index]).Im() * (bus[indexer].V[temp_index]).Im();
02735 tempIcalcImag += (bus[indexer].full_Y_load[temp_index]).Re() * (bus[indexer].V[temp_index]).Im() + (bus[indexer].full_Y_load[temp_index]).Im() * (bus[indexer].V[temp_index]).Re();
02736 }
02737
02738
02739
02740 temp_index_b = -1;
02741 switch(bus[indexer].phases & 0x07) {
02742 case 0x01:
02743 {
02744 temp_index_b=2;
02745 break;
02746 }
02747 case 0x02:
02748 {
02749 temp_index_b=1;
02750 break;
02751 }
02752 case 0x03:
02753 {
02754 if (jindex==0)
02755 temp_index_b=1;
02756 else
02757 temp_index_b=2;
02758 break;
02759 }
02760 case 0x04:
02761 {
02762 temp_index_b=0;
02763 break;
02764 }
02765 case 0x05:
02766 {
02767 if (jindex==0)
02768 temp_index_b=0;
02769 else
02770 temp_index_b=2;
02771 break;
02772 }
02773 case 0x06:
02774 {
02775 if (jindex==0)
02776 temp_index_b=0;
02777 else
02778 temp_index_b=1;
02779 break;
02780 }
02781 case 0x07:
02782 {
02783 temp_index_b = jindex;
02784 break;
02785 }
02786 }
02787
02788 if (temp_index_b==-1)
02789 {
02790 GL_THROW("NR: A voltage index failed to be found.");
02791 }
02792
02793 for (kindexer=0; kindexer<(bus[indexer].Link_Table_Size); kindexer++)
02794 {
02795
02796 jindexer=bus[indexer].Link_Table[kindexer];
02797
02798 if (branch[jindexer].from == indexer)
02799 {
02800
02801 if ((branch[jindexer].phases & 0x80) == 0x80)
02802 {
02803 proceed_flag = false;
02804 phase_worka = branch[jindexer].phases & 0x07;
02805
02806 if (kindex==0)
02807 {
02808 switch (bus[indexer].phases & 0x07) {
02809 case 0x01:
02810 {
02811 if (phase_worka==0x01)
02812 proceed_flag=true;
02813 break;
02814 }
02815 case 0x02:
02816 {
02817 if (phase_worka==0x02)
02818 proceed_flag=true;
02819 break;
02820 }
02821 case 0x03:
02822 {
02823 if ((jindex==0) && (phase_worka==0x02))
02824 proceed_flag=true;
02825 else if ((jindex==1) && (phase_worka==0x01))
02826 proceed_flag=true;
02827 else
02828 ;
02829 break;
02830 }
02831 case 0x04:
02832 {
02833 if (phase_worka==0x04)
02834 proceed_flag=true;
02835 break;
02836 }
02837 case 0x05:
02838 {
02839 if ((jindex==0) && (phase_worka==0x04))
02840 proceed_flag=true;
02841 else if ((jindex==1) && (phase_worka==0x01))
02842 proceed_flag=true;
02843 else
02844 ;
02845 break;
02846 }
02847 case 0x06:
02848 case 0x07:
02849 {
02850 if ((jindex==0) && (phase_worka==0x04))
02851 proceed_flag=true;
02852 else if ((jindex==1) && (phase_worka==0x02))
02853 proceed_flag=true;
02854 else if ((jindex==2) && (phase_worka==0x01))
02855 proceed_flag=true;
02856 else;
02857 break;
02858 }
02859 }
02860 }
02861
02862 if (proceed_flag)
02863 {
02864 work_vals_char_0 = temp_index_b*3;
02865
02866
02867 tempIcalcReal += (-(branch[jindexer].Yfrom[work_vals_char_0])).Re() * (bus[branch[jindexer].to].V[0]).Re() - (-(branch[jindexer].Yfrom[work_vals_char_0])).Im() * (bus[branch[jindexer].to].V[0]).Im();
02868 tempIcalcImag += (-(branch[jindexer].Yfrom[work_vals_char_0])).Re() * (bus[branch[jindexer].to].V[0]).Im() + (-(branch[jindexer].Yfrom[work_vals_char_0])).Im() * (bus[branch[jindexer].to].V[0]).Re();
02869
02870
02871 tempIcalcReal += (-(branch[jindexer].Yfrom[work_vals_char_0+1])).Re() * (bus[branch[jindexer].to].V[1]).Re() - (-(branch[jindexer].Yfrom[work_vals_char_0+1])).Im() * (bus[branch[jindexer].to].V[1]).Im();
02872 tempIcalcImag += (-(branch[jindexer].Yfrom[work_vals_char_0+1])).Re() * (bus[branch[jindexer].to].V[1]).Im() + (-(branch[jindexer].Yfrom[work_vals_char_0+1])).Im() * (bus[branch[jindexer].to].V[1]).Re();
02873
02874 }
02875 }
02876 else
02877 {
02878 work_vals_char_0 = temp_index_b*3+temp_index;
02879 work_vals_double_0 = (-branch[jindexer].Yfrom[work_vals_char_0]).Re();
02880 work_vals_double_1 = (-branch[jindexer].Yfrom[work_vals_char_0]).Im();
02881 work_vals_double_2 = (bus[branch[jindexer].to].V[temp_index]).Re();
02882 work_vals_double_3 = (bus[branch[jindexer].to].V[temp_index]).Im();
02883
02884 tempIcalcReal += work_vals_double_0 * work_vals_double_2 - work_vals_double_1 * work_vals_double_3;
02885 tempIcalcImag += work_vals_double_0 * work_vals_double_3 + work_vals_double_1 * work_vals_double_2;
02886
02887 }
02888 }
02889 if (branch[jindexer].to == indexer)
02890 {
02891 work_vals_char_0 = temp_index_b*3+temp_index;
02892 work_vals_double_0 = (-branch[jindexer].Yto[work_vals_char_0]).Re();
02893 work_vals_double_1 = (-branch[jindexer].Yto[work_vals_char_0]).Im();
02894 work_vals_double_2 = (bus[branch[jindexer].from].V[temp_index]).Re();
02895 work_vals_double_3 = (bus[branch[jindexer].from].V[temp_index]).Im();
02896
02897 tempIcalcReal += work_vals_double_0 * work_vals_double_2 - work_vals_double_1 * work_vals_double_3;
02898 tempIcalcImag += work_vals_double_0 * work_vals_double_3 + work_vals_double_1 * work_vals_double_2;
02899 }
02900 else;
02901 }
02902 }
02903
02904
02905 if ((bus[indexer].type > 1) && (bus[indexer].swing_functions_enabled == true))
02906 {
02907 if ((powerflow_type == PF_DYNINIT) && (bus[indexer].DynCurrent != NULL))
02908 {
02909
02910 if (bus[indexer].full_Y != NULL)
02911 {
02912
02913 temp_complex_2 = bus[indexer].V[jindex] * complex(tempIcalcReal,-tempIcalcImag);
02914
02915 if (powerflow_values->island_matrix_values[island_loop_index].iteration_count>0)
02916 {
02917
02918 (*bus[indexer].PGenTotal) += temp_complex_2 - temp_complex_3/3.0;
02919 }
02920
02921
02922 work_vals_double_0 = (bus[indexer].V[temp_index_b]).Mag()*(bus[indexer].V[temp_index_b]).Mag();
02923
02924 if (work_vals_double_0!=0)
02925 {
02926 work_vals_double_1 = (bus[indexer].V[temp_index_b]).Re();
02927 work_vals_double_2 = (bus[indexer].V[temp_index_b]).Im();
02928 work_vals_double_3 = (tempPbus * work_vals_double_1 + tempQbus * work_vals_double_2)/ (work_vals_double_0) - tempIcalcReal;
02929 work_vals_double_4 = (tempPbus * work_vals_double_2 - tempQbus * work_vals_double_1)/ (work_vals_double_0) - tempIcalcImag;
02930
02931
02932 work_vals_double_0 = sqrt((work_vals_double_3*work_vals_double_3+work_vals_double_4*work_vals_double_4));
02933 }
02934 else
02935 {
02936 work_vals_double_0 = 0.0;
02937 }
02938
02939 if (work_vals_double_0 > bus[indexer].max_volt_error)
02940 {
02941 powerflow_values->island_matrix_values[island_loop_index].swing_converged=false;
02942 }
02943 }
02944 else
02945 {
02946
02947 work_vals_double_0 = (bus[indexer].V[temp_index_b]).Mag()*(bus[indexer].V[temp_index_b]).Mag();
02948
02949 if (work_vals_double_0!=0)
02950 {
02951 work_vals_double_1 = (bus[indexer].V[temp_index_b]).Re();
02952 work_vals_double_2 = (bus[indexer].V[temp_index_b]).Im();
02953 work_vals_double_3 = (tempPbus * work_vals_double_1 + tempQbus * work_vals_double_2)/ (work_vals_double_0) - tempIcalcReal;
02954 work_vals_double_4 = (tempPbus * work_vals_double_2 - tempQbus * work_vals_double_1)/ (work_vals_double_0) - tempIcalcImag;
02955
02956
02957 work_vals_double_0 = sqrt((work_vals_double_3*work_vals_double_3+work_vals_double_4*work_vals_double_4));
02958 }
02959 else
02960 {
02961 work_vals_double_0 = 0.0;
02962 }
02963
02964 if (work_vals_double_0 > bus[indexer].max_volt_error)
02965 {
02966 powerflow_values->island_matrix_values[island_loop_index].swing_converged=false;
02967 }
02968
02969
02970 bus[indexer].DynCurrent[jindex] = complex(tempIcalcReal,tempIcalcImag);
02971 }
02972 }
02973
02974
02975
02976 if (NR_busdata[indexer].BusHistTerm != NULL)
02977 {
02978 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + jindex] = NR_busdata[indexer].BusHistTerm[jindex].Re();
02979 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = NR_busdata[indexer].BusHistTerm[jindex].Im();
02980 }
02981 else
02982 {
02983 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + jindex] = 0.0;
02984 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = 0.0;
02985 }
02986
02987
02988 }
02989 else
02990 {
02991 work_vals_double_0 = (bus[indexer].V[temp_index_b]).Mag()*(bus[indexer].V[temp_index_b]).Mag();
02992
02993 if (work_vals_double_0!=0)
02994 {
02995 work_vals_double_1 = (bus[indexer].V[temp_index_b]).Re();
02996 work_vals_double_2 = (bus[indexer].V[temp_index_b]).Im();
02997
02998
02999 if (NR_busdata[indexer].BusHistTerm != NULL)
03000 {
03001 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+ powerflow_values->BA_diag[indexer].size + jindex] = (tempPbus * work_vals_double_1 + tempQbus * work_vals_double_2)/ (work_vals_double_0) + NR_busdata[indexer].BusHistTerm[jindex].Re() - tempIcalcReal ;
03002 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = (tempPbus * work_vals_double_2 - tempQbus * work_vals_double_1)/ (work_vals_double_0) + NR_busdata[indexer].BusHistTerm[jindex].Im() - tempIcalcImag;
03003 }
03004 else
03005 {
03006 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+ powerflow_values->BA_diag[indexer].size + jindex] = (tempPbus * work_vals_double_1 + tempQbus * work_vals_double_2)/ (work_vals_double_0) - tempIcalcReal ;
03007 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = (tempPbus * work_vals_double_2 - tempQbus * work_vals_double_1)/ (work_vals_double_0) - tempIcalcImag;
03008 }
03009
03010
03011 if (NR_busdata[indexer].BusSatTerm != NULL)
03012 {
03013 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+ powerflow_values->BA_diag[indexer].size + jindex] -= NR_busdata[indexer].BusSatTerm[jindex].Re();
03014 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] -= NR_busdata[indexer].BusSatTerm[jindex].Im();
03015 }
03016 }
03017 else
03018 {
03019 if (NR_busdata[indexer].BusHistTerm != NULL)
03020 {
03021 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + jindex] = NR_busdata[indexer].BusHistTerm[jindex].Re();
03022 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = NR_busdata[indexer].BusHistTerm[jindex].Im();
03023 }
03024 else
03025 {
03026 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + jindex] = 0.0;
03027 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] = 0.0;
03028 }
03029
03030
03031 if (NR_busdata[indexer].BusSatTerm != NULL)
03032 {
03033 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+ powerflow_values->BA_diag[indexer].size + jindex] -= NR_busdata[indexer].BusSatTerm[jindex].Re();
03034 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + jindex] -= NR_busdata[indexer].BusSatTerm[jindex].Im();
03035 }
03036 }
03037 }
03038 }
03039 }
03040
03041
03042 if (powerflow_type != PF_NORMAL)
03043 {
03044
03045 if ((*bus[indexer].dynamics_enabled==true) && (bus[indexer].DynCurrent != NULL))
03046 {
03047
03048 if (bus[indexer].full_Y != NULL)
03049 {
03050 if ((bus[indexer].phases & 0x07) == 0x07)
03051 {
03052
03053 if ((powerflow_values->island_matrix_values[island_loop_index].swing_is_a_swing == true) && (powerflow_type == PF_DYNINIT))
03054 {
03055
03056 temp_complex_0 += complex((*bus[indexer].PGenTotal).Re(),-(*bus[indexer].PGenTotal).Im());
03057
03058
03059 if (temp_complex_1.Mag() > 0.0)
03060 {
03061 temp_complex_2 = temp_complex_0/temp_complex_1;
03062
03063
03064 bus[indexer].DynCurrent[0] = temp_complex_2;
03065 bus[indexer].DynCurrent[1] = temp_complex_2*avalsq;
03066 bus[indexer].DynCurrent[2] = temp_complex_2*aval;
03067 }
03068 else
03069 {
03070 bus[indexer].DynCurrent[0] = complex(0.0,0.0);
03071 bus[indexer].DynCurrent[1] = complex(0.0,0.0);
03072 bus[indexer].DynCurrent[2] = complex(0.0,0.0);
03073 }
03074 }
03075
03076
03077 if ((bus[indexer].type == 0) || ((bus[indexer].type>1) && (bus[indexer].swing_functions_enabled == false)))
03078 {
03079
03080 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size] += bus[indexer].DynCurrent[0].Re();
03081 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + 1] += bus[indexer].DynCurrent[1].Re();
03082 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + 2] += bus[indexer].DynCurrent[2].Re();
03083
03084 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc] += bus[indexer].DynCurrent[0].Im();
03085 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + 1] += bus[indexer].DynCurrent[1].Im();
03086 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + 2] += bus[indexer].DynCurrent[2].Im();
03087 }
03088
03089 }
03090 else
03091 {
03092
03093 if ((bus[indexer].origphases & 0x07) == 0x07)
03094 {
03095
03096 bus[indexer].DynCurrent[0] = complex(0.0,0.0);
03097 bus[indexer].DynCurrent[1] = complex(0.0,0.0);
03098 bus[indexer].DynCurrent[2] = complex(0.0,0.0);
03099 }
03100 else
03101 {
03102 GL_THROW("NR_solver: bus:%s has tried updating deltamode dynamics, but is not three-phase!",bus[indexer].name);
03103
03104
03105
03106
03107
03108 }
03109 }
03110 }
03111 else
03112 {
03113
03114
03115 if ((bus[indexer].type == 0) || ((bus[indexer].type>1) && (bus[indexer].swing_functions_enabled == false)))
03116 {
03117
03118 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size] += bus[indexer].DynCurrent[0].Re();
03119 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + 1] += bus[indexer].DynCurrent[1].Re();
03120 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size + 2] += bus[indexer].DynCurrent[2].Re();
03121
03122 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc] += bus[indexer].DynCurrent[0].Im();
03123 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + 1] += bus[indexer].DynCurrent[1].Im();
03124 powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[2*bus[indexer].Matrix_Loc + 2] += bus[indexer].DynCurrent[2].Im();
03125 }
03126
03127 }
03128 }
03129 }
03130 }
03131
03132 }
03133
03134
03135 compute_load_values(bus_count,bus,powerflow_values,true,island_loop_index);
03136
03137
03138
03139
03140 powerflow_values->island_matrix_values[island_loop_index].size_diag_update = 0;
03141
03142
03143 for (jindexer=0; jindexer<bus_count;jindexer++)
03144 {
03145
03146 if (bus[jindexer].island_number == island_loop_index)
03147 {
03148 if (bus[jindexer].type != 1)
03149 powerflow_values->island_matrix_values[island_loop_index].size_diag_update += powerflow_values->BA_diag[jindexer].size;
03150
03151 }
03152
03153 }
03154
03155 if (powerflow_values->island_matrix_values[island_loop_index].Y_diag_update == NULL)
03156 {
03157 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update = (Y_NR *)gl_malloc((4*powerflow_values->island_matrix_values[island_loop_index].size_diag_update) *sizeof(Y_NR));
03158
03159
03160 if (powerflow_values->island_matrix_values[island_loop_index].Y_diag_update == NULL)
03161 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
03162
03163
03164 powerflow_values->island_matrix_values[island_loop_index].max_size_diag_update = powerflow_values->island_matrix_values[island_loop_index].size_diag_update;
03165 }
03166 else if (powerflow_values->island_matrix_values[island_loop_index].size_diag_update > powerflow_values->island_matrix_values[island_loop_index].max_size_diag_update)
03167 {
03168
03169 gl_free(powerflow_values->island_matrix_values[island_loop_index].Y_diag_update);
03170
03171
03172 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update = (Y_NR *)gl_malloc((4*powerflow_values->island_matrix_values[island_loop_index].size_diag_update) *sizeof(Y_NR));
03173
03174
03175 if (powerflow_values->island_matrix_values[island_loop_index].Y_diag_update == NULL)
03176 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
03177
03178
03179 powerflow_values->island_matrix_values[island_loop_index].max_size_diag_update = powerflow_values->island_matrix_values[island_loop_index].size_diag_update;
03180
03181
03182 powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed = true;
03183 }
03184
03185 indexer = 0;
03186
03187 for (jindexer=0; jindexer<bus_count; jindexer++)
03188 {
03189
03190 if (bus[jindexer].island_number == island_loop_index)
03191 {
03192 if ((bus[jindexer].type > 1) && (bus[jindexer].swing_functions_enabled == true))
03193 {
03194 for (jindex=0; jindex<powerflow_values->BA_diag[jindexer].size; jindex++)
03195 {
03196 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex;
03197 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind;
03198 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = 1e10;
03199 indexer += 1;
03200
03201 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex;
03202 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind + powerflow_values->BA_diag[jindexer].size;
03203 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][jindex]).Re();
03204 indexer += 1;
03205
03206 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex + powerflow_values->BA_diag[jindexer].size;
03207 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind - powerflow_values->BA_diag[jindexer].size;
03208 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][jindex]).Re();
03209 indexer += 1;
03210
03211 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex + powerflow_values->BA_diag[jindexer].size;
03212 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind;
03213 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = 1e10;
03214 indexer += 1;
03215 }
03216 }
03217
03218 if ((bus[jindexer].type == 0) || ((bus[jindexer].type > 1) && bus[jindexer].swing_functions_enabled == false))
03219 {
03220 for (jindex=0; jindex<powerflow_values->BA_diag[jindexer].size; jindex++)
03221 {
03222 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex;
03223 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind;
03224 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][jindex]).Im() + bus[jindexer].Jacob_A[jindex];
03225 indexer += 1;
03226
03227 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex;
03228 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind + powerflow_values->BA_diag[jindexer].size;
03229 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][jindex]).Re() + bus[jindexer].Jacob_B[jindex];
03230 indexer += 1;
03231
03232 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex + powerflow_values->BA_diag[jindexer].size;
03233 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = 2*bus[jindexer].Matrix_Loc + jindex;
03234 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = (powerflow_values->BA_diag[jindexer].Y[jindex][jindex]).Re() + bus[jindexer].Jacob_C[jindex];
03235 indexer += 1;
03236
03237 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind = 2*bus[jindexer].Matrix_Loc + jindex + powerflow_values->BA_diag[jindexer].size;
03238 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind;
03239 powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value = -(powerflow_values->BA_diag[jindexer].Y[jindex][jindex]).Im() + bus[jindexer].Jacob_D[jindex];
03240 indexer += 1;
03241 }
03242 }
03243 }
03244 }
03245
03246
03247 powerflow_values->island_matrix_values[island_loop_index].size_Amatrix = powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ*2 + powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed*2 + 4*powerflow_values->island_matrix_values[island_loop_index].size_diag_update;
03248
03249
03250 if (powerflow_values->island_matrix_values[island_loop_index].size_Amatrix==0)
03251 {
03252
03253 if (NR_islands_detected > 1)
03254 {
03255 gl_warning("Empty powerflow connectivity matrix for island %d, your system is empty!",(island_loop_index+1));
03256
03257
03258
03259
03260
03261 }
03262 else
03263 {
03264 gl_warning("Empty powerflow connectivity matrix, your system is empty!");
03265
03266
03267
03268
03269
03270 }
03271
03272
03273 powerflow_values->island_matrix_values[island_loop_index].return_code = 0;
03274
03275
03276
03277
03278 if (island_loop_index >= (NR_islands_detected - 1))
03279 {
03280
03281 still_iterating_islands = false;
03282 }
03283 else
03284 {
03285
03286 island_loop_index++;
03287
03288
03289 still_iterating_islands = true;
03290 }
03291
03292
03293 continue;
03294 }
03295
03296 if (powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix == NULL)
03297 {
03298 powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix = (SPARSE*) gl_malloc(sizeof(SPARSE));
03299
03300
03301 if (powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix == NULL)
03302 GL_THROW("NR: Failed to allocate memory for one of the necessary matrices");
03303
03304
03305 sparse_init(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix, powerflow_values->island_matrix_values[island_loop_index].size_Amatrix, 6*powerflow_values->island_matrix_values[island_loop_index].bus_count);
03306 }
03307 else if (powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed)
03308 {
03309
03310 sparse_clear(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix);
03311
03312
03313 sparse_init(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix, powerflow_values->island_matrix_values[island_loop_index].size_Amatrix, 6*powerflow_values->island_matrix_values[island_loop_index].bus_count);
03314 }
03315 else
03316 {
03317
03318 sparse_reset(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix, 6*powerflow_values->island_matrix_values[island_loop_index].bus_count);
03319 }
03320
03321
03322 for (indexer=0; indexer<powerflow_values->island_matrix_values[island_loop_index].size_offdiag_PQ*2; indexer++)
03323 {
03324 row = powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ[indexer].row_ind;
03325 col = powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ[indexer].col_ind;
03326 value = powerflow_values->island_matrix_values[island_loop_index].Y_offdiag_PQ[indexer].Y_value;
03327 sparse_add(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix, row, col, value, bus, bus_count, powerflow_values, island_loop_index);
03328 }
03329
03330
03331 for (indexer=0; indexer< (powerflow_values->island_matrix_values[island_loop_index].size_diag_fixed*2); indexer++)
03332 {
03333 row = powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed[indexer].row_ind;
03334 col = powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed[indexer].col_ind;
03335 value = powerflow_values->island_matrix_values[island_loop_index].Y_diag_fixed[indexer].Y_value;
03336 sparse_add(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix, row, col, value, bus, bus_count, powerflow_values, island_loop_index);
03337 }
03338
03339
03340 for (indexer=0; indexer< (4*powerflow_values->island_matrix_values[island_loop_index].size_diag_update); indexer++)
03341 {
03342 row = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].row_ind;
03343 col = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].col_ind;
03344 value = powerflow_values->island_matrix_values[island_loop_index].Y_diag_update[indexer].Y_value;
03345 sparse_add(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix, row, col, value, bus, bus_count, powerflow_values, island_loop_index);
03346 }
03347
03348
03349 if (NRMatDumpMethod != MD_NONE)
03350 {
03351
03352
03353
03354 something_has_been_output = false;
03355
03356
03357 if ((NRMatDumpMethod == MD_ALL) || ((NRMatDumpMethod != MD_ALL) && (powerflow_values->island_matrix_values[island_loop_index].iteration_count == 0)))
03358 {
03359
03360 FPoutVal=fopen(MDFileName,"at");
03361
03362
03363 if ((NRMatReferences == true) && (powerflow_values->island_matrix_values[island_loop_index].iteration_count == 0))
03364 {
03365
03366 if (NR_islands_detected > 1)
03367 {
03368 fprintf(FPoutVal,"Matrix Index information for this call in island:%d - start,stop,name\n",(island_loop_index+1));
03369 }
03370 else
03371 {
03372 fprintf(FPoutVal,"Matrix Index information for this call - start,stop,name\n");
03373 }
03374
03375 for (indexer=0; indexer<bus_count; indexer++)
03376 {
03377
03378 if (bus[indexer].island_number == island_loop_index)
03379 {
03380
03381 jindexer = 2*bus[indexer].Matrix_Loc;
03382 kindexer = jindexer + 2*powerflow_values->BA_diag[indexer].size - 1;
03383
03384
03385 fprintf(FPoutVal,"%d,%d,%s\n",jindexer,kindexer,bus[indexer].name);
03386
03387
03388 something_has_been_output = true;
03389 }
03390 }
03391
03392
03393 fprintf(FPoutVal,"\n");
03394 }
03395
03396
03397 fprintf(FPoutVal,"Timestamp: %lld - Iteration %lld\n",gl_globalclock,powerflow_values->island_matrix_values[island_loop_index].iteration_count);
03398
03399
03400 if (something_has_been_output == true)
03401 {
03402
03403 fprintf(FPoutVal,"Matrix Information - non-zero element count = %d\n",powerflow_values->island_matrix_values[island_loop_index].size_Amatrix);
03404
03405
03406
03407
03408 fprintf(FPoutVal,"Matrix Information - row, column, value\n");
03409
03410
03411 temp_element = NULL;
03412
03413
03414 for (jindexer=0; jindexer<powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix->ncols; jindexer++)
03415 {
03416
03417 temp_element = powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix->cols[jindexer];
03418
03419
03420 if (temp_element != NULL)
03421 {
03422
03423 fprintf(FPoutVal,"%d,%d,%f\n",temp_element->row_ind,jindexer,temp_element->value);
03424
03425
03426 while (temp_element->next != NULL)
03427 {
03428
03429 temp_element = temp_element->next;
03430
03431
03432 fprintf(FPoutVal,"%d,%d,%f\n",temp_element->row_ind,jindexer,temp_element->value);
03433 }
03434 }
03435
03436 }
03437
03438
03439 fprintf(FPoutVal,"\n");
03440 }
03441 else
03442 {
03443
03444 fprintf(FPoutVal,"Empty island (possibly powerflow-removed)\n\n");
03445 }
03446
03447
03448 fclose(FPoutVal);
03449
03450
03451 if (NRMatDumpMethod == MD_ONCE)
03452 {
03453
03454 if (NR_islands_detected > 1)
03455 {
03456
03457 if (island_loop_index >= (NR_islands_detected - 1))
03458 {
03459 NRMatDumpMethod = MD_NONE;
03460 }
03461
03462 }
03463 else
03464 {
03465 NRMatDumpMethod = MD_NONE;
03466 }
03467 }
03468 }
03469 }
03470
03472 m = 2*powerflow_values->island_matrix_values[island_loop_index].total_variables;
03473 n = 2*powerflow_values->island_matrix_values[island_loop_index].total_variables;
03474 nnz = powerflow_values->island_matrix_values[island_loop_index].size_Amatrix;
03475
03476 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.a_LU == NULL)
03477 {
03478
03479 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.a_LU = (double *) gl_malloc(nnz *sizeof(double));
03480 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.a_LU==NULL)
03481 {
03482 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03483
03484
03485
03486
03487
03488 }
03489
03490 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rows_LU = (int *) gl_malloc(nnz *sizeof(int));
03491 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rows_LU == NULL)
03492 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03493
03494 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.cols_LU = (int *) gl_malloc((n+1) *sizeof(int));
03495 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.cols_LU == NULL)
03496 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03497
03498
03499 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU = (double *) gl_malloc(m *sizeof(double));
03500 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU == NULL)
03501 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03502
03503 if (matrix_solver_method==MM_SUPERLU)
03504 {
03506 curr_island_superLU_vars->perm_r = (int *) gl_malloc(m *sizeof(int));
03507 if (curr_island_superLU_vars->perm_r == NULL)
03508 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03509
03510 curr_island_superLU_vars->perm_c = (int *) gl_malloc(n *sizeof(int));
03511 if (curr_island_superLU_vars->perm_c == NULL)
03512 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03513
03514
03515 curr_island_superLU_vars->A_LU.Store = (void *)gl_malloc(sizeof(NCformat));
03516 if (curr_island_superLU_vars->A_LU.Store == NULL)
03517 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03518
03519 curr_island_superLU_vars->B_LU.Store = (void *)gl_malloc(sizeof(DNformat));
03520 if (curr_island_superLU_vars->B_LU.Store == NULL)
03521 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03522
03523
03524 curr_island_superLU_vars->A_LU.Stype = SLU_NC;
03525 curr_island_superLU_vars->A_LU.Dtype = SLU_D;
03526 curr_island_superLU_vars->A_LU.Mtype = SLU_GE;
03527 curr_island_superLU_vars->A_LU.nrow = n;
03528 curr_island_superLU_vars->A_LU.ncol = m;
03529
03530
03531 curr_island_superLU_vars->B_LU.Stype = SLU_DN;
03532 curr_island_superLU_vars->B_LU.Dtype = SLU_D;
03533 curr_island_superLU_vars->B_LU.Mtype = SLU_GE;
03534 curr_island_superLU_vars->B_LU.nrow = m;
03535 curr_island_superLU_vars->B_LU.ncol = 1;
03536 }
03537 else if (matrix_solver_method == MM_EXTERN)
03538 {
03539
03540 ((void (*)(void *,unsigned int, unsigned int, bool))(LUSolverFcns.ext_alloc))(powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars,n,n,NR_admit_change);
03541 }
03542 else
03543 {
03544 GL_THROW("Invalid matrix solution method specified for NR solver!");
03545
03546 }
03547
03548
03549 powerflow_values->island_matrix_values[island_loop_index].prev_m = m;
03550 }
03551 else if (powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed)
03552 {
03553
03554 gl_free(powerflow_values->island_matrix_values[island_loop_index].matrices_LU.a_LU);
03555 gl_free(powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rows_LU);
03556 gl_free(powerflow_values->island_matrix_values[island_loop_index].matrices_LU.cols_LU);
03557 gl_free(powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU);
03558
03559 if (matrix_solver_method==MM_SUPERLU)
03560 {
03561
03562 gl_free(curr_island_superLU_vars->perm_r);
03563 gl_free(curr_island_superLU_vars->perm_c);
03564 }
03565
03566
03567
03568 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.a_LU = (double *) gl_malloc(nnz *sizeof(double));
03569 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.a_LU==NULL)
03570 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03571
03572 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rows_LU = (int *) gl_malloc(nnz *sizeof(int));
03573 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rows_LU == NULL)
03574 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03575
03576 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.cols_LU = (int *) gl_malloc((n+1) *sizeof(int));
03577 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.cols_LU == NULL)
03578 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03579
03580
03581 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU = (double *) gl_malloc(m *sizeof(double));
03582 if (powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU == NULL)
03583 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03584
03585 if (matrix_solver_method==MM_SUPERLU)
03586 {
03588 curr_island_superLU_vars->perm_r = (int *) gl_malloc(m *sizeof(int));
03589 if (curr_island_superLU_vars->perm_r == NULL)
03590 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03591
03592 curr_island_superLU_vars->perm_c = (int *) gl_malloc(n *sizeof(int));
03593 if (curr_island_superLU_vars->perm_c == NULL)
03594 GL_THROW("NR: One of the SuperLU solver matrices failed to allocate");
03595
03596
03597 curr_island_superLU_vars->A_LU.Stype = SLU_NC;
03598 curr_island_superLU_vars->A_LU.Dtype = SLU_D;
03599 curr_island_superLU_vars->A_LU.Mtype = SLU_GE;
03600 curr_island_superLU_vars->A_LU.nrow = n;
03601 curr_island_superLU_vars->A_LU.ncol = m;
03602
03603
03604 curr_island_superLU_vars->B_LU.Stype = SLU_DN;
03605 curr_island_superLU_vars->B_LU.Dtype = SLU_D;
03606 curr_island_superLU_vars->B_LU.Mtype = SLU_GE;
03607 curr_island_superLU_vars->B_LU.nrow = m;
03608 curr_island_superLU_vars->B_LU.ncol = 1;
03609 }
03610 else if (matrix_solver_method == MM_EXTERN)
03611 {
03612
03613 ((void (*)(void *,unsigned int, unsigned int, bool))(LUSolverFcns.ext_alloc))(powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars,n,n,NR_admit_change);
03614 }
03615 else
03616 {
03617 GL_THROW("Invalid matrix solution method specified for NR solver!");
03618
03619 }
03620
03621
03622 powerflow_values->island_matrix_values[island_loop_index].prev_m = m;
03623 }
03624 else if (powerflow_values->island_matrix_values[island_loop_index].prev_m != m)
03625 {
03626 if (matrix_solver_method==MM_SUPERLU)
03627 {
03628
03629 curr_island_superLU_vars->A_LU.nrow = n;
03630 curr_island_superLU_vars->A_LU.ncol = m;
03631
03632 curr_island_superLU_vars->B_LU.nrow = m;
03633 }
03634 else if (matrix_solver_method == MM_EXTERN)
03635 {
03636
03637 ((void (*)(void *,unsigned int, unsigned int, bool))(LUSolverFcns.ext_alloc))(powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars,n,n,NR_admit_change);
03638 }
03639 else
03640 {
03641 GL_THROW("Invalid matrix solution method specified for NR solver!");
03642
03643 }
03644
03645
03646 powerflow_values->island_matrix_values[island_loop_index].prev_m = m;
03647 }
03648
03649 #ifndef MT
03650 if (matrix_solver_method==MM_SUPERLU)
03651 {
03652
03653 set_default_options ( &options );
03654 }
03655
03656 #endif
03657
03658 sparse_tonr(powerflow_values->island_matrix_values[island_loop_index].Y_Amatrix, &powerflow_values->island_matrix_values[island_loop_index].matrices_LU);
03659 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.cols_LU[n] = nnz ;
03660
03661
03662 if (mesh_imped_vals == NULL)
03663 {
03664 for (temp_index_c=0;temp_index_c<m;temp_index_c++)
03665 {
03666 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU[temp_index_c] = powerflow_values->island_matrix_values[island_loop_index].deltaI_NR[temp_index_c];
03667 }
03668 }
03669
03670
03671 if (matrix_solver_method==MM_SUPERLU)
03672 {
03674
03675 Astore = (NCformat*)curr_island_superLU_vars->A_LU.Store;
03676 Astore->nnz = nnz;
03677 Astore->nzval = powerflow_values->island_matrix_values[island_loop_index].matrices_LU.a_LU;
03678 Astore->rowind = powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rows_LU;
03679 Astore->colptr = powerflow_values->island_matrix_values[island_loop_index].matrices_LU.cols_LU;
03680
03681
03682
03683 Bstore = (DNformat*)curr_island_superLU_vars->B_LU.Store;
03684 Bstore->lda = m;
03685 Bstore->nzval = powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU;
03686
03687
03688
03689 if (mesh_imped_vals != NULL)
03690 {
03691
03692 mesh_imped_vals->return_code = 0;
03693
03694
03695 if (mesh_imped_vals->NodeRefNum > 0)
03696 {
03697 tempa = 2*bus[mesh_imped_vals->NodeRefNum].Matrix_Loc;
03698 }
03699 else
03700 {
03701 gl_error("solver_nr: Node reference for mesh fault calculation was invalid -- it may be a child");
03702
03703
03704
03705
03706
03707
03708
03709 mesh_imped_vals->return_code = 0;
03710 *bad_computations = true;
03711
03712
03713 NR_solver_working = false;
03714
03715
03716 return 0;
03717 }
03718
03719
03720 temp_size = 2*powerflow_values->BA_diag[mesh_imped_vals->NodeRefNum].size;
03721
03722
03723 if (temp_size == 0)
03724 {
03725 gl_error("solver_nr: Mesh fault impedance extract failed - invalid node");
03726
03727
03728
03729
03730
03731
03732
03733
03734 *bad_computations = true;
03735 mesh_imped_vals->return_code = 0;
03736
03737
03738 NR_solver_working = false;
03739
03740
03741 return 0;
03742 }
03743
03744
03745
03746 for (jindex=0; jindex<9; jindex++)
03747 {
03748 mesh_imped_vals->z_matrix[jindex] = complex(0.0,0.0);
03749 }
03750
03751
03752 for (jindex=0; jindex<6; jindex++)
03753 {
03754 for (kindex=0; kindex<6; kindex++)
03755 {
03756 temp_z_store[jindex][kindex] = 0.0;
03757 }
03758 }
03759
03760
03761 for (kindex=0; kindex<temp_size; kindex++)
03762 {
03763
03764
03765 for (temp_index_c=0;temp_index_c<m;temp_index_c++)
03766 {
03767 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU[temp_index_c] = 0.0;
03768 }
03769
03770
03771 powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU[tempa + kindex] = 1.0;
03772
03773
03774 #ifdef MT
03775
03776
03777
03778 get_perm_c(1, &curr_island_superLU_vars->A_LU, curr_island_superLU_vars->perm_c);
03779
03780
03781 pdgssv(NR_superLU_procs, &curr_island_superLU_vars->A_LU, curr_island_superLU_vars->perm_c, curr_island_superLU_vars->perm_r, &L_LU, &U_LU, &curr_island_superLU_vars->B_LU, &info);
03782
03783
03784
03785 Destroy_SuperNode_SCP(&L_LU);
03786 Destroy_CompCol_NCP(&U_LU);
03787 #else
03788
03789
03790 StatInit ( &stat );
03791
03792
03793 dgssv(&options, &curr_island_superLU_vars->A_LU, curr_island_superLU_vars->perm_c, curr_island_superLU_vars->perm_r, &L_LU, &U_LU, &curr_island_superLU_vars->B_LU, &stat, &info);
03794
03795
03796
03797 Destroy_SuperNode_Matrix( &L_LU );
03798 Destroy_CompCol_Matrix( &U_LU );
03799 StatFree ( &stat );
03800 #endif
03801
03802 if (info != 0)
03803 {
03804 gl_error("solver_nr: superLU failed mesh fault matrix inversion with code %d",info);
03805
03806
03807
03808
03809
03810
03811
03812 *bad_computations = true;
03813 mesh_imped_vals->return_code = 0;
03814
03815
03816 NR_solver_working = false;
03817
03818
03819 return 0;
03820 }
03821
03822
03823
03824 sol_LU = (double*) ((DNformat*) curr_island_superLU_vars->B_LU.Store)->nzval;
03825
03826
03827 for (jindex=0; jindex<temp_size; jindex++)
03828 {
03829 temp_z_store[jindex][kindex] = sol_LU[tempa+jindex];
03830 }
03831 }
03832
03833
03834 switch(bus[mesh_imped_vals->NodeRefNum].phases & 0x07) {
03835 case 0x01:
03836 {
03837 mesh_imped_vals->z_matrix[8] = complex(temp_z_store[1][0],temp_z_store[1][1]);
03838 break;
03839 }
03840 case 0x02:
03841 {
03842 mesh_imped_vals->z_matrix[4] = complex(temp_z_store[1][0],temp_z_store[1][1]);
03843 break;
03844 }
03845 case 0x03:
03846 {
03847
03848 mesh_imped_vals->z_matrix[4] = complex(temp_z_store[2][0],temp_z_store[2][2]);
03849 mesh_imped_vals->z_matrix[5] = complex(temp_z_store[2][1],temp_z_store[2][3]);
03850
03851
03852 mesh_imped_vals->z_matrix[7] = complex(temp_z_store[3][0],temp_z_store[3][2]);
03853 mesh_imped_vals->z_matrix[8] = complex(temp_z_store[3][1],temp_z_store[3][3]);
03854 break;
03855 }
03856 case 0x04:
03857 {
03858 mesh_imped_vals->z_matrix[0] = complex(temp_z_store[1][0],temp_z_store[1][1]);
03859 break;
03860 }
03861 case 0x05:
03862 {
03863
03864 mesh_imped_vals->z_matrix[0] = complex(temp_z_store[2][0],temp_z_store[2][2]);
03865 mesh_imped_vals->z_matrix[2] = complex(temp_z_store[2][1],temp_z_store[2][3]);
03866
03867
03868 mesh_imped_vals->z_matrix[6] = complex(temp_z_store[3][0],temp_z_store[3][2]);
03869 mesh_imped_vals->z_matrix[8] = complex(temp_z_store[3][1],temp_z_store[3][3]);
03870 break;
03871 }
03872 case 0x06:
03873 {
03874
03875 mesh_imped_vals->z_matrix[0] = complex(temp_z_store[2][0],temp_z_store[2][2]);
03876 mesh_imped_vals->z_matrix[1] = complex(temp_z_store[2][1],temp_z_store[2][3]);
03877
03878
03879 mesh_imped_vals->z_matrix[3] = complex(temp_z_store[3][0],temp_z_store[3][2]);
03880 mesh_imped_vals->z_matrix[4] = complex(temp_z_store[3][1],temp_z_store[3][3]);
03881 break;
03882 }
03883 case 0x07:
03884 {
03885
03886 mesh_imped_vals->z_matrix[0] = complex(temp_z_store[3][0],temp_z_store[3][3]);
03887 mesh_imped_vals->z_matrix[1] = complex(temp_z_store[3][1],temp_z_store[3][4]);
03888 mesh_imped_vals->z_matrix[2] = complex(temp_z_store[3][2],temp_z_store[3][5]);
03889
03890
03891 mesh_imped_vals->z_matrix[3] = complex(temp_z_store[4][0],temp_z_store[4][3]);
03892 mesh_imped_vals->z_matrix[4] = complex(temp_z_store[4][1],temp_z_store[4][4]);
03893 mesh_imped_vals->z_matrix[5] = complex(temp_z_store[4][2],temp_z_store[4][5]);
03894
03895
03896 mesh_imped_vals->z_matrix[6] = complex(temp_z_store[5][0],temp_z_store[5][3]);
03897 mesh_imped_vals->z_matrix[7] = complex(temp_z_store[5][1],temp_z_store[5][4]);
03898 mesh_imped_vals->z_matrix[8] = complex(temp_z_store[5][2],temp_z_store[5][5]);
03899 break;
03900 }
03901 default:
03902 {
03903 gl_error("solver_nr: Mesh fault impedance extract failed - invalid node");
03904
03905
03906
03907 *bad_computations = true;
03908 mesh_imped_vals->return_code = 0;
03909
03910
03911 NR_solver_working = false;
03912
03913
03914 return 0;
03915 }
03916 }
03917
03918
03919 *bad_computations = false;
03920 mesh_imped_vals->return_code = 1;
03921
03922
03923 return 1;
03924 }
03925 else
03926 {
03927 #ifdef MT
03928
03929
03930
03931 get_perm_c(1, &curr_island_superLU_vars->A_LU, curr_island_superLU_vars->perm_c);
03932
03933
03934 pdgssv(NR_superLU_procs, &curr_island_superLU_vars->A_LU, curr_island_superLU_vars->perm_c, curr_island_superLU_vars->perm_r, &L_LU, &U_LU, &curr_island_superLU_vars->B_LU, &powerflow_values->island_matrix_values[island_loop_index].solver_info);
03935 #else
03936
03937
03938 StatInit ( &stat );
03939
03940
03941 dgssv(&options, &curr_island_superLU_vars->A_LU, curr_island_superLU_vars->perm_c, curr_island_superLU_vars->perm_r, &L_LU, &U_LU, &curr_island_superLU_vars->B_LU, &stat, &powerflow_values->island_matrix_values[island_loop_index].solver_info);
03942 #endif
03943
03944 sol_LU = (double*) ((DNformat*) curr_island_superLU_vars->B_LU.Store)->nzval;
03945 }
03946 }
03947 else if (matrix_solver_method==MM_EXTERN)
03948 {
03949
03950 if (mesh_imped_vals != NULL)
03951 {
03952 gl_error("solver_nr: Mesh impedance attempted from unsupported LU solver");
03953
03954
03955
03956
03957
03958
03959
03960 mesh_imped_vals->return_code = 2;
03961
03962
03963 ((void (*)(void *, bool))(LUSolverFcns.ext_destroy))(powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars,powerflow_values->island_matrix_values[island_loop_index].new_iteration_required);
03964
03965
03966 *bad_computations = true;
03967
03968
03969 NR_solver_working = false;
03970
03971
03972 return 0;
03973 }
03974
03975
03976
03977 powerflow_values->island_matrix_values[island_loop_index].solver_info = ((int (*)(void *,NR_SOLVER_VARS *, unsigned int, unsigned int))(LUSolverFcns.ext_solve))(powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars,&powerflow_values->island_matrix_values[island_loop_index].matrices_LU,n,1);
03978
03979
03980 sol_LU = powerflow_values->island_matrix_values[island_loop_index].matrices_LU.rhs_LU;
03981 }
03982 else
03983 {
03984 GL_THROW("Invalid matrix solution method specified for NR solver!");
03985
03986
03987
03988
03989 }
03990
03991
03992 powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge = 0;
03993
03994 temp_index = -1;
03995 temp_index_b = -1;
03996 powerflow_values->island_matrix_values[island_loop_index].new_iteration_required = false;
03997
03998 for (indexer=0; indexer<bus_count; indexer++)
03999 {
04000
04001 if (bus[indexer].island_number == island_loop_index)
04002 {
04003
04004 if ((bus[indexer].type == 0) || ((bus[indexer].type > 1) && (bus[indexer].swing_functions_enabled == false)))
04005 {
04006
04007 if ((bus[indexer].phases & 0x80) == 0x80)
04008 {
04009
04010 DVConvCheck[0]=complex(sol_LU[2*bus[indexer].Matrix_Loc],sol_LU[(2*bus[indexer].Matrix_Loc+2)]);
04011 DVConvCheck[1]=complex(sol_LU[(2*bus[indexer].Matrix_Loc+1)],sol_LU[(2*bus[indexer].Matrix_Loc+3)]);
04012 bus[indexer].V[0] += DVConvCheck[0];
04013 bus[indexer].V[1] += DVConvCheck[1];
04014
04015
04016 CurrConvVal=DVConvCheck[0].Mag();
04017 if (CurrConvVal > powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge)
04018 powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge=CurrConvVal;
04019
04020 if (CurrConvVal > bus[indexer].max_volt_error)
04021 powerflow_values->island_matrix_values[island_loop_index].new_iteration_required=true;
04022
04023 CurrConvVal=DVConvCheck[1].Mag();
04024 if (CurrConvVal > powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge)
04025 powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge=CurrConvVal;
04026
04027 if (CurrConvVal > bus[indexer].max_volt_error)
04028 powerflow_values->island_matrix_values[island_loop_index].new_iteration_required=true;
04029 }
04030 else
04031 {
04032 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
04033 {
04034 switch(bus[indexer].phases & 0x07) {
04035 case 0x01:
04036 {
04037 temp_index=0;
04038 temp_index_b=2;
04039 break;
04040 }
04041 case 0x02:
04042 {
04043 temp_index=0;
04044 temp_index_b=1;
04045 break;
04046 }
04047 case 0x03:
04048 {
04049 if (jindex==0)
04050 {
04051 temp_index=0;
04052 temp_index_b=1;
04053 }
04054 else
04055 {
04056 temp_index=1;
04057 temp_index_b=2;
04058 }
04059
04060 break;
04061 }
04062 case 0x04:
04063 {
04064 temp_index=0;
04065 temp_index_b=0;
04066 break;
04067 }
04068 case 0x05:
04069 {
04070 if (jindex==0)
04071 {
04072 temp_index=0;
04073 temp_index_b=0;
04074 }
04075 else
04076 {
04077 temp_index=1;
04078 temp_index_b=2;
04079 }
04080 break;
04081 }
04082 case 0x06:
04083 case 0x07:
04084 {
04085 temp_index = jindex;
04086 temp_index_b = jindex;
04087 break;
04088 }
04089 }
04090
04091 if ((temp_index==-1) || (temp_index_b==-1))
04092 {
04093 GL_THROW("NR: An error occurred indexing voltage updates");
04094
04095
04096
04097
04098
04099 }
04100
04101 DVConvCheck[jindex]=complex(sol_LU[(2*bus[indexer].Matrix_Loc+temp_index)],sol_LU[(2*bus[indexer].Matrix_Loc+powerflow_values->BA_diag[indexer].size+temp_index)]);
04102 bus[indexer].V[temp_index_b] += DVConvCheck[jindex];
04103
04104
04105 CurrConvVal=DVConvCheck[jindex].Mag();
04106 if (CurrConvVal > bus[indexer].max_volt_error)
04107 powerflow_values->island_matrix_values[island_loop_index].new_iteration_required=true;
04108
04109 if (CurrConvVal > powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge)
04110 powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge = CurrConvVal;
04111
04112 }
04113 }
04114 }
04115
04116
04117
04118 if ((bus[indexer].ExtraCurrentInjFunc != NULL) && ((bus[indexer].type == 0) || ((bus[indexer].type > 1) && (bus[indexer].swing_functions_enabled == false))))
04119 {
04120
04121 call_return_status = ((STATUS (*)(OBJECT *))(*bus[indexer].ExtraCurrentInjFunc))(bus[indexer].ExtraCurrentInjFuncObject);
04122
04123
04124 if (call_return_status == FAILED)
04125 {
04126 GL_THROW("External current injection update failed for device %s",bus[indexer].ExtraCurrentInjFuncObject->name ? bus[indexer].ExtraCurrentInjFuncObject->name : "Unnamed");
04127
04128
04129
04130
04131 }
04132
04133 }
04134 }
04135
04136 }
04137
04138
04139
04140 if ((enable_inrush_calculations == true) && (deltatimestep_running > 0))
04141 {
04142
04143 powerflow_values->island_matrix_values[island_loop_index].SaturationMismatchPresent = false;
04144
04145
04146 for (indexer=0; indexer<branch_count; indexer++)
04147 {
04148
04149 if (branch[indexer].island_number == island_loop_index)
04150 {
04151
04152 if (branch[indexer].ExtraDeltaModeFunc != NULL)
04153 {
04154
04155 func_result_val = ((int (*)(OBJECT *,bool *))(*branch[indexer].ExtraDeltaModeFunc))(branch[indexer].obj,&(powerflow_values->island_matrix_values[island_loop_index].SaturationMismatchPresent));
04156
04157
04158 if (func_result_val != 1)
04159 {
04160 GL_THROW("Extra delta update failed for device %s",branch[indexer].name ? branch[indexer].name : "Unnamed");
04161
04162
04163
04164
04165 }
04166
04167 }
04168
04169 }
04170
04171 }
04172 }
04173
04174
04175 if (powerflow_type == PF_DYNINIT)
04176 {
04177
04178 if (powerflow_values->island_matrix_values[island_loop_index].new_iteration_required == false)
04179 {
04180 if (powerflow_values->island_matrix_values[island_loop_index].swing_converged == false)
04181 {
04182
04183 if (NR_islands_detected > 1)
04184 {
04185
04186 gl_verbose("NR_Solver: swing bus failed balancing criteria on island %d - making it a PQ for future dynamics",(island_loop_index+1));
04187 }
04188 else
04189 {
04190
04191 gl_verbose("NR_Solver: swing bus failed balancing criteria - making it a PQ for future dynamics");
04192 }
04193
04194
04195 if (powerflow_values->island_matrix_values[island_loop_index].swing_is_a_swing == true)
04196 {
04197
04198 powerflow_values->island_matrix_values[island_loop_index].swing_is_a_swing = false;
04199
04200
04201 for (indexer=0; indexer<bus_count; indexer++)
04202 {
04203
04204 if (bus[indexer].island_number == island_loop_index)
04205 {
04206
04207 if ((bus[indexer].type > 1) && (bus[indexer].swing_functions_enabled == true))
04208 {
04209
04210 if ((*bus[indexer].dynamics_enabled==true) && (bus[indexer].DynCurrent != NULL))
04211 {
04212
04213 bus[indexer].swing_functions_enabled = false;
04214
04215
04216 powerflow_values->island_matrix_values[island_loop_index].new_iteration_required=true;
04217 }
04218 }
04219
04220 }
04221
04222 }
04223 }
04224
04225 }
04226
04227 }
04228
04229 }
04230
04231
04232 if (powerflow_values->island_matrix_values[island_loop_index].SaturationMismatchPresent == true)
04233 {
04234
04235 powerflow_values->island_matrix_values[island_loop_index].new_iteration_required = true;
04236 }
04237
04238
04239
04240 powerflow_values->island_matrix_values[island_loop_index].NR_realloc_needed = false;
04241
04242 if (matrix_solver_method==MM_SUPERLU)
04243 {
04244
04245 #ifdef MT
04246
04247 Destroy_SuperNode_SCP(&L_LU);
04248 Destroy_CompCol_NCP(&U_LU);
04249 #else
04250
04251 Destroy_SuperNode_Matrix( &L_LU );
04252 Destroy_CompCol_Matrix( &U_LU );
04253 StatFree ( &stat );
04254 #endif
04255 }
04256 else if (matrix_solver_method==MM_EXTERN)
04257 {
04258
04259 ((void (*)(void *, bool))(LUSolverFcns.ext_destroy))(powerflow_values->island_matrix_values[island_loop_index].LU_solver_vars,powerflow_values->island_matrix_values[island_loop_index].new_iteration_required);
04260 }
04261 else
04262 {
04263 GL_THROW("Invalid matrix solution method specified for NR solver!");
04264
04265 }
04266
04267
04268 if (powerflow_values->island_matrix_values[island_loop_index].solver_info != 0)
04269 {
04270 powerflow_values->island_matrix_values[island_loop_index].new_iteration_required = false;
04271 }
04272
04273
04274 proceed_to_next_island = false;
04275
04276
04277 if (powerflow_values->island_matrix_values[island_loop_index].new_iteration_required == true)
04278 {
04279
04280 if (powerflow_values->island_matrix_values[island_loop_index].iteration_count>=NR_iteration_limit)
04281 {
04282 powerflow_values->island_matrix_values[island_loop_index].return_code = -(powerflow_values->island_matrix_values[island_loop_index].iteration_count);
04283
04284 if (NR_islands_detected > 1)
04285 {
04286 gl_verbose("Max solver mismatch of failed solution %f on island %d\n",powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge,(island_loop_index+1));
04287 }
04288 else
04289 {
04290 gl_verbose("Max solver mismatch of failed solution %f\n",powerflow_values->island_matrix_values[island_loop_index].max_mismatch_converge);
04291 }
04292
04293
04294 proceed_to_next_island = true;
04295 }
04296 else
04297 {
04298 powerflow_values->island_matrix_values[island_loop_index].iteration_count++;
04299
04300
04301 proceed_to_next_island = false;
04302 }
04303 }
04304 else
04305 {
04306
04307 proceed_to_next_island = true;
04308
04309
04310 if (powerflow_values->island_matrix_values[island_loop_index].solver_info == 0)
04311 {
04312
04313 if (NR_islands_detected > 1)
04314 {
04315 gl_verbose("Power flow calculation for island %d converges at Iteration %d \n",(island_loop_index+1),(powerflow_values->island_matrix_values[island_loop_index].iteration_count+1));
04316 }
04317 else
04318 {
04319 gl_verbose("Power flow calculation converges at Iteration %d \n",(powerflow_values->island_matrix_values[island_loop_index].iteration_count+1));
04320 }
04321
04322
04323 powerflow_values->island_matrix_values[island_loop_index].return_code = powerflow_values->island_matrix_values[island_loop_index].iteration_count;
04324 }
04325 else
04326 {
04327
04328
04329
04330 if (NR_islands_detected > 1)
04331 {
04332
04333 if (matrix_solver_method==MM_SUPERLU)
04334 {
04335 gl_verbose("superLU failed out of island %d with return value %d",(island_loop_index+1),powerflow_values->island_matrix_values[island_loop_index].solver_info);
04336 }
04337 else if (matrix_solver_method==MM_EXTERN)
04338 {
04339 gl_verbose("External LU solver failed out of island %d with return value %d",(island_loop_index+1),powerflow_values->island_matrix_values[island_loop_index].solver_info);
04340 }
04341
04342 }
04343 else
04344 {
04345
04346 if (matrix_solver_method==MM_SUPERLU)
04347 {
04348 gl_verbose("superLU failed out with return value %d",powerflow_values->island_matrix_values[island_loop_index].solver_info);
04349 }
04350 else if (matrix_solver_method==MM_EXTERN)
04351 {
04352 gl_verbose("External LU solver failed out with return value %d",powerflow_values->island_matrix_values[island_loop_index].solver_info);
04353 }
04354
04355 }
04356
04357
04358 powerflow_values->island_matrix_values[island_loop_index].return_code = 0;
04359 }
04360 }
04361
04362
04363 if (proceed_to_next_island == true)
04364 {
04365
04366 if (island_loop_index >= (NR_islands_detected - 1))
04367 {
04368
04369 still_iterating_islands = false;
04370 }
04371 else
04372 {
04373
04374 island_loop_index++;
04375
04376
04377 still_iterating_islands = true;
04378 }
04379 }
04380 }
04381
04382
04383
04384 *bad_computations = true;
04385 return_value_for_solver_NR = 0;
04386
04387
04388 for (island_loop_index=0; island_loop_index<NR_islands_detected; island_loop_index++)
04389 {
04390 if (powerflow_values->island_matrix_values[island_loop_index].return_code >= 0)
04391 {
04392
04393 if (powerflow_values->island_matrix_values[island_loop_index].new_iteration_required == false)
04394 {
04395 *bad_computations = false;
04396
04397
04398 if (return_value_for_solver_NR >= 0)
04399 {
04400 if (powerflow_values->island_matrix_values[island_loop_index].iteration_count > return_value_for_solver_NR)
04401 {
04402 return_value_for_solver_NR = powerflow_values->island_matrix_values[island_loop_index].iteration_count;
04403 }
04404 }
04405
04406 }
04407 else if (NR_island_fail_method == true)
04408 {
04409 *bad_computations = false;
04410
04411
04412
04413 return_value_for_solver_NR = -NR_iteration_limit;
04414
04415
04416 temp_fxn_val = (FUNCTIONADDR)(gl_get_function(fault_check_object,"island_removal_function"));
04417
04418
04419 if (temp_fxn_val == NULL)
04420 {
04421 GL_THROW("NR: Unable to map island removal function");
04422
04423
04424
04425
04426
04427 }
04428
04429
04430 call_return_status = ((STATUS (*)(OBJECT *,int))(temp_fxn_val))(fault_check_object,island_loop_index);
04431
04432
04433 if (call_return_status != SUCCESS)
04434 {
04435 GL_THROW("NR: Failed to remove island %d from the powerflow",(island_loop_index+1));
04436
04437
04438
04439
04440 }
04441
04442 }
04443
04444 }
04445 else
04446 {
04447 *bad_computations = false;
04448
04449
04450 if (NR_island_fail_method == true)
04451 {
04452
04453
04454 return_value_for_solver_NR = -NR_iteration_limit;
04455
04456
04457 temp_fxn_val = (FUNCTIONADDR)(gl_get_function(fault_check_object,"island_removal_function"));
04458
04459
04460 if (temp_fxn_val == NULL)
04461 {
04462 GL_THROW("NR: Unable to map island removal function");
04463
04464 }
04465
04466
04467 call_return_status = ((STATUS (*)(OBJECT *,int))(temp_fxn_val))(fault_check_object,island_loop_index);
04468
04469
04470 if (call_return_status != SUCCESS)
04471 {
04472 GL_THROW("NR: Failed to remove island %d from the powerflow",(island_loop_index+1));
04473
04474 }
04475
04476 }
04477 else
04478 {
04479
04480 if (powerflow_values->island_matrix_values[island_loop_index].return_code < return_value_for_solver_NR)
04481 {
04482 return_value_for_solver_NR = powerflow_values->island_matrix_values[island_loop_index].return_code;
04483 }
04484 }
04485 }
04486
04487 }
04488
04489
04490 NR_solver_working = false;
04491
04492
04493 if (*bad_computations == true)
04494 {
04495
04496 return 0;
04497 }
04498 else
04499 {
04500
04501 return return_value_for_solver_NR;
04502 }
04503 }
04504
04505
04506
04507
04508
04509 void compute_load_values(unsigned int bus_count, BUSDATA *bus, NR_SOLVER_STRUCT *powerflow_values, bool jacobian_pass, int island_number)
04510 {
04511 unsigned int indexer;
04512 double adjust_nominal_voltage_val, adjust_nominal_voltaged_val;
04513 double tempPbus, tempQbus;
04514 double adjust_temp_voltage_mag[6];
04515 complex adjust_temp_nominal_voltage[6], adjusted_constant_current[6];
04516 complex delta_current[3], voltageDel[3], undeltacurr[3];
04517 complex temp_current[3], temp_store[3];
04518 char jindex, temp_index, temp_index_b;
04519
04520
04521 for (indexer=0; indexer<bus_count; indexer++)
04522 {
04523
04524 if (bus[indexer].island_number == island_number)
04525 {
04526 if ((bus[indexer].phases & 0x08) == 0x08)
04527 {
04528
04529 if (*bus[indexer].dynamics_enabled == true)
04530 {
04531
04532 adjust_nominal_voltage_val = bus[indexer].volt_base * sqrt(3.0);
04533
04534
04535 adjust_temp_nominal_voltage[0].SetPolar(adjust_nominal_voltage_val,PI/6.0);
04536 adjust_temp_nominal_voltage[1].SetPolar(adjust_nominal_voltage_val,-1.0*PI/2.0);
04537 adjust_temp_nominal_voltage[2].SetPolar(adjust_nominal_voltage_val,5.0*PI/6.0);
04538
04539
04540 voltageDel[0] = bus[indexer].V[0] - bus[indexer].V[1];
04541 voltageDel[1] = bus[indexer].V[1] - bus[indexer].V[2];
04542 voltageDel[2] = bus[indexer].V[2] - bus[indexer].V[0];
04543
04544
04545 adjust_temp_voltage_mag[0] = voltageDel[0].Mag();
04546 adjust_temp_voltage_mag[1] = voltageDel[1].Mag();
04547 adjust_temp_voltage_mag[2] = voltageDel[2].Mag();
04548
04549
04550 if ((bus[indexer].I[0] != 0.0) && (adjust_temp_voltage_mag[0] != 0.0))
04551 {
04552
04553 adjusted_constant_current[0] = ~(adjust_temp_nominal_voltage[0] * ~bus[indexer].I[0] * adjust_temp_voltage_mag[0] / (voltageDel[0] * adjust_nominal_voltage_val));
04554 }
04555 else
04556 {
04557 adjusted_constant_current[0] = complex(0.0,0.0);
04558 }
04559
04560
04561 if ((bus[indexer].I[1] != 0.0) && (adjust_temp_voltage_mag[1] != 0.0))
04562 {
04563
04564 adjusted_constant_current[1] = ~(adjust_temp_nominal_voltage[1] * ~bus[indexer].I[1] * adjust_temp_voltage_mag[1] / (voltageDel[1] * adjust_nominal_voltage_val));
04565 }
04566 else
04567 {
04568 adjusted_constant_current[1] = complex(0.0,0.0);
04569 }
04570
04571
04572 if ((bus[indexer].I[2] != 0.0) && (adjust_temp_voltage_mag[2] != 0.0))
04573 {
04574
04575 adjusted_constant_current[2] = ~(adjust_temp_nominal_voltage[2] * ~bus[indexer].I[2] * adjust_temp_voltage_mag[2] / (voltageDel[2] * adjust_nominal_voltage_val));
04576 }
04577 else
04578 {
04579 adjusted_constant_current[2] = complex(0.0,0.0);
04580 }
04581
04582
04583 if ((bus[indexer].phases & 0x10) == 0x10)
04584 {
04585
04586 adjust_nominal_voltage_val = bus[indexer].volt_base;
04587
04588
04589 adjust_temp_nominal_voltage[3].SetPolar(bus[indexer].volt_base,0.0);
04590 adjust_temp_nominal_voltage[4].SetPolar(bus[indexer].volt_base,-2.0*PI/3.0);
04591 adjust_temp_nominal_voltage[5].SetPolar(bus[indexer].volt_base,2.0*PI/3.0);
04592
04593
04594 adjust_temp_voltage_mag[3] = bus[indexer].V[0].Mag();
04595 adjust_temp_voltage_mag[4] = bus[indexer].V[1].Mag();
04596 adjust_temp_voltage_mag[5] = bus[indexer].V[2].Mag();
04597
04598
04599 if ((bus[indexer].extra_var[6] != 0.0) && (adjust_temp_voltage_mag[3] != 0.0))
04600 {
04601
04602 adjusted_constant_current[3] = ~(adjust_temp_nominal_voltage[3] * ~bus[indexer].extra_var[6] * adjust_temp_voltage_mag[3] / (bus[indexer].V[0] * adjust_nominal_voltage_val));
04603 }
04604 else
04605 {
04606 adjusted_constant_current[3] = complex(0.0,0.0);
04607 }
04608
04609
04610 if ((bus[indexer].extra_var[7] != 0.0) && (adjust_temp_voltage_mag[4] != 0.0))
04611 {
04612
04613 adjusted_constant_current[4] = ~(adjust_temp_nominal_voltage[4] * ~bus[indexer].extra_var[7] * adjust_temp_voltage_mag[4] / (bus[indexer].V[1] * adjust_nominal_voltage_val));
04614 }
04615 else
04616 {
04617 adjusted_constant_current[4] = complex(0.0,0.0);
04618 }
04619
04620
04621 if ((bus[indexer].extra_var[8] != 0.0) && (adjust_temp_voltage_mag[5] != 0.0))
04622 {
04623
04624 adjusted_constant_current[5] = ~(adjust_temp_nominal_voltage[5] * ~bus[indexer].extra_var[8] * adjust_temp_voltage_mag[5] / (bus[indexer].V[2] * adjust_nominal_voltage_val));
04625 }
04626 else
04627 {
04628 adjusted_constant_current[5] = complex(0.0,0.0);
04629 }
04630 }
04631 else
04632 {
04633
04634 adjusted_constant_current[3] = complex(0.0,0.0);
04635 adjusted_constant_current[4] = complex(0.0,0.0);
04636 adjusted_constant_current[5] = complex(0.0,0.0);
04637 }
04638 }
04639 else
04640 {
04641 adjusted_constant_current[0] = bus[indexer].I[0];
04642 adjusted_constant_current[1] = bus[indexer].I[1];
04643 adjusted_constant_current[2] = bus[indexer].I[2];
04644
04645
04646 if ((bus[indexer].phases & 0x10) == 0x10)
04647 {
04648
04649 adjusted_constant_current[3] = bus[indexer].extra_var[6];
04650 adjusted_constant_current[4] = bus[indexer].extra_var[7];
04651 adjusted_constant_current[5] = bus[indexer].extra_var[8];
04652 }
04653 else
04654 {
04655 adjusted_constant_current[3] = complex(0.0,0.0);
04656 adjusted_constant_current[4] = complex(0.0,0.0);
04657 adjusted_constant_current[5] = complex(0.0,0.0);
04658 }
04659 }
04660
04661
04662 if ((bus[indexer].phases & 0x06) == 0x06)
04663 {
04664
04665 voltageDel[0] = bus[indexer].V[0] - bus[indexer].V[1];
04666
04667
04668 delta_current[0] = (voltageDel[0] == 0) ? 0 : ~(bus[indexer].S[0]/voltageDel[0]);
04669
04670
04671 delta_current[0] += voltageDel[0] * (bus[indexer].Y[0]);
04672 }
04673 else
04674 {
04675
04676 voltageDel[0] = complex(0.0,0.0);
04677 delta_current[0] = complex(0.0,0.0);
04678 }
04679
04680 if ((bus[indexer].phases & 0x03) == 0x03)
04681 {
04682
04683 voltageDel[1] = bus[indexer].V[1] - bus[indexer].V[2];
04684
04685
04686 delta_current[1] = (voltageDel[1] == 0) ? 0 : ~(bus[indexer].S[1]/voltageDel[1]);
04687
04688
04689 delta_current[1] += voltageDel[1] * (bus[indexer].Y[1]);
04690 }
04691 else
04692 {
04693
04694 voltageDel[1] = complex(0.0,0.0);
04695 delta_current[1] = complex(0.0,0.0);
04696 }
04697
04698 if ((bus[indexer].phases & 0x05) == 0x05)
04699 {
04700
04701 voltageDel[2] = bus[indexer].V[2] - bus[indexer].V[0];
04702
04703
04704 delta_current[2] = (voltageDel[2] == 0) ? 0 : ~(bus[indexer].S[2]/voltageDel[2]);
04705
04706
04707 delta_current[2] += voltageDel[2] * (bus[indexer].Y[2]);
04708 }
04709 else
04710 {
04711
04712 voltageDel[2] = complex(0.0,0.0);
04713 delta_current[2] = complex(0.0,0.0);
04714 }
04715
04716
04717
04718 if ((bus[indexer].phases & 0x04) == 0x04)
04719 {
04720 undeltacurr[0]=(adjusted_constant_current[0]+delta_current[0])-(adjusted_constant_current[2]+delta_current[2]);
04721
04722
04723 if ((bus[indexer].phases & 0x10) == 0x10)
04724 {
04725
04726 undeltacurr[0] += (bus[indexer].V[0] == 0) ? 0 : ~(bus[indexer].extra_var[0]/bus[indexer].V[0]);
04727
04728
04729 undeltacurr[0] += bus[indexer].extra_var[3]*bus[indexer].V[0];
04730
04731
04732 undeltacurr[0] += adjusted_constant_current[3];
04733 }
04734 }
04735 else
04736 {
04737
04738 undeltacurr[0] = complex(0.0,0.0);
04739 }
04740
04741 if ((bus[indexer].phases & 0x02) == 0x02)
04742 {
04743 undeltacurr[1]=(adjusted_constant_current[1]+delta_current[1])-(adjusted_constant_current[0]+delta_current[0]);
04744
04745
04746 if ((bus[indexer].phases & 0x10) == 0x10)
04747 {
04748
04749 undeltacurr[1] += (bus[indexer].V[1] == 0) ? 0 : ~(bus[indexer].extra_var[1]/bus[indexer].V[1]);
04750
04751
04752 undeltacurr[1] += bus[indexer].extra_var[4]*bus[indexer].V[1];
04753
04754
04755 undeltacurr[1] += adjusted_constant_current[4];
04756 }
04757 }
04758 else
04759 {
04760
04761 undeltacurr[1] = complex(0.0,0.0);
04762 }
04763
04764 if ((bus[indexer].phases & 0x01) == 0x01)
04765 {
04766 undeltacurr[2]=(adjusted_constant_current[2]+delta_current[2])-(adjusted_constant_current[1]+delta_current[1]);
04767
04768
04769 if ((bus[indexer].phases & 0x10) == 0x10)
04770 {
04771
04772 undeltacurr[2] += (bus[indexer].V[2] == 0) ? 0 : ~(bus[indexer].extra_var[2]/bus[indexer].V[2]);
04773
04774
04775 undeltacurr[2] += bus[indexer].extra_var[5]*bus[indexer].V[2];
04776
04777
04778 undeltacurr[2] += adjusted_constant_current[5];
04779 }
04780 }
04781 else
04782 {
04783
04784 undeltacurr[2] = complex(0.0,0.0);
04785 }
04786
04787
04788
04789 temp_index = -1;
04790 temp_index_b = -1;
04791
04792 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
04793 {
04794 switch(bus[indexer].phases & 0x07) {
04795 case 0x01:
04796 {
04797 temp_index=0;
04798 temp_index_b=2;
04799 break;
04800 }
04801 case 0x02:
04802 {
04803 temp_index=0;
04804 temp_index_b=1;
04805 break;
04806 }
04807 case 0x03:
04808 {
04809 if (jindex==0)
04810 {
04811 temp_index=0;
04812 temp_index_b=1;
04813 }
04814 else
04815 {
04816 temp_index=1;
04817 temp_index_b=2;
04818 }
04819 break;
04820 }
04821 case 0x04:
04822 {
04823 temp_index=0;
04824 temp_index_b=0;
04825 break;
04826 }
04827 case 0x05:
04828 {
04829 if (jindex==0)
04830 {
04831 temp_index=0;
04832 temp_index_b=0;
04833 }
04834 else
04835 {
04836 temp_index=1;
04837 temp_index_b=2;
04838 }
04839 break;
04840 }
04841 case 0x06:
04842 case 0x07:
04843 {
04844 temp_index=jindex;
04845 temp_index_b=jindex;
04846 break;
04847 }
04848 default:
04849 break;
04850 }
04851
04852 if (jacobian_pass == false)
04853 {
04854 if ((temp_index==-1) || (temp_index_b==-1))
04855 {
04856 GL_THROW("NR: A scheduled power update element failed.");
04857
04858 }
04859
04860
04861 tempPbus = (undeltacurr[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Re() + (undeltacurr[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Im();
04862 bus[indexer].PL[temp_index] = tempPbus;
04863
04864
04865 tempQbus = (undeltacurr[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Im() - (undeltacurr[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Re();
04866 bus[indexer].QL[temp_index] = tempQbus;
04867 }
04868 else
04869 {
04870 if ((temp_index==-1) || (temp_index_b==-1))
04871 {
04872 GL_THROW("NR: A Jacobian update element failed.");
04873
04874 }
04875
04876 if ((bus[indexer].V[temp_index_b]).Mag()!=0)
04877 {
04878 bus[indexer].Jacob_A[temp_index] = ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Re() + (undeltacurr[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
04879 bus[indexer].Jacob_B[temp_index] = -((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Im() + (undeltacurr[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
04880 bus[indexer].Jacob_C[temp_index] =((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Im() - (undeltacurr[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
04881 bus[indexer].Jacob_D[temp_index] = ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Re() - (undeltacurr[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
04882 }
04883 else
04884 {
04885 bus[indexer].Jacob_A[temp_index] = -1e-4;
04886 bus[indexer].Jacob_B[temp_index] = -1e-4;
04887 bus[indexer].Jacob_C[temp_index] = -1e-4;
04888 bus[indexer].Jacob_D[temp_index] = -1e-4;
04889 }
04890 }
04891 }
04892 }
04893 else if ((bus[indexer].phases & 0x80) == 0x80)
04894 {
04895
04896
04897 voltageDel[0] = bus[indexer].V[0] + bus[indexer].V[1];
04898
04899
04900 temp_current[0] = bus[indexer].I[0];
04901 temp_current[1] = bus[indexer].I[1];
04902 temp_current[2] = *bus[indexer].extra_var;
04903
04904
04905
04906 if ((bus[indexer].prerot_I[2] != 0.0) && (*bus[indexer].dynamics_enabled == true))
04907 temp_current[2] += bus[indexer].prerot_I[2];
04908
04909
04910 temp_current[0] += bus[indexer].V[0] == 0.0 ? 0.0 : ~(bus[indexer].S[0]/bus[indexer].V[0]);
04911 temp_current[1] += bus[indexer].V[1] == 0.0 ? 0.0 : ~(bus[indexer].S[1]/bus[indexer].V[1]);
04912 temp_current[2] += voltageDel[0] == 0.0 ? 0.0 : ~(bus[indexer].S[2]/voltageDel[0]);
04913
04914
04915 temp_current[0] += bus[indexer].Y[0]*bus[indexer].V[0];
04916 temp_current[1] += bus[indexer].Y[1]*bus[indexer].V[1];
04917 temp_current[2] += bus[indexer].Y[2]*voltageDel[0];
04918
04919
04920 if ((bus[indexer].phases & 0x40) == 0x40)
04921 {
04922
04923 temp_store[0].SetPolar(1.0,bus[indexer].V[0].Arg());
04924 temp_store[1].SetPolar(1.0,bus[indexer].V[1].Arg());
04925 temp_store[2].SetPolar(1.0,voltageDel[0].Arg());
04926
04927
04928 delta_current[0] = bus[indexer].house_var[0]/(~temp_store[0]);
04929 delta_current[1] = bus[indexer].house_var[1]/(~temp_store[1]);
04930 delta_current[2] = bus[indexer].house_var[2]/(~temp_store[2]);
04931
04932
04933 temp_current[0] += delta_current[0];
04934 temp_current[1] += delta_current[1];
04935 temp_current[2] += delta_current[2];
04936 }
04937
04938 if (jacobian_pass == false)
04939 {
04940
04941 temp_store[0] = temp_current[0] + temp_current[2];
04942 temp_store[1] = -temp_current[1] - temp_current[2];
04943
04944
04945 bus[indexer].PL[0] = temp_store[0].Re();
04946 bus[indexer].QL[0] = temp_store[0].Im();
04947
04948 bus[indexer].PL[1] = temp_store[1].Re();
04949 bus[indexer].QL[1] = temp_store[1].Im();
04950 }
04951 else
04952 {
04953
04954 temp_store[0] = -(temp_current[0] + temp_current[2]);
04955 temp_store[1] = -(-temp_current[1] - temp_current[2]);
04956
04957 for (jindex=0; jindex<2; jindex++)
04958 {
04959 if ((bus[indexer].V[jindex]).Mag()!=0)
04960 {
04961 bus[indexer].Jacob_A[jindex] = ((bus[indexer].V[jindex]).Re()*(bus[indexer].V[jindex]).Im()*(temp_store[jindex]).Re() + (temp_store[jindex]).Im() *pow((bus[indexer].V[jindex]).Im(),2))/pow((bus[indexer].V[jindex]).Mag(),3);
04962 bus[indexer].Jacob_B[jindex] = -((bus[indexer].V[jindex]).Re()*(bus[indexer].V[jindex]).Im()*(temp_store[jindex]).Im() + (temp_store[jindex]).Re() *pow((bus[indexer].V[jindex]).Re(),2))/pow((bus[indexer].V[jindex]).Mag(),3);
04963 bus[indexer].Jacob_C[jindex] =((bus[indexer].V[jindex]).Re()*(bus[indexer].V[jindex]).Im()*(temp_store[jindex]).Im() - (temp_store[jindex]).Re() *pow((bus[indexer].V[jindex]).Im(),2))/pow((bus[indexer].V[jindex]).Mag(),3);
04964 bus[indexer].Jacob_D[jindex] = ((bus[indexer].V[jindex]).Re()*(bus[indexer].V[jindex]).Im()*(temp_store[jindex]).Re() - (temp_store[jindex]).Im() *pow((bus[indexer].V[jindex]).Re(),2))/pow((bus[indexer].V[jindex]).Mag(),3);
04965 }
04966 else
04967 {
04968 bus[indexer].Jacob_A[jindex]= -1e-4;
04969 bus[indexer].Jacob_B[jindex]= -1e-4;
04970 bus[indexer].Jacob_C[jindex]= -1e-4;
04971 bus[indexer].Jacob_D[jindex]= -1e-4;
04972 }
04973 }
04974
04975
04976 bus[indexer].Jacob_A[2] = 0.0;
04977 bus[indexer].Jacob_B[2] = 0.0;
04978 bus[indexer].Jacob_C[2] = 0.0;
04979 bus[indexer].Jacob_D[2] = 0.0;
04980 }
04981 }
04982 else
04983 {
04984
04985 if (*bus[indexer].dynamics_enabled == true)
04986 {
04987
04988 adjust_nominal_voltage_val = bus[indexer].volt_base;
04989
04990
04991 adjust_temp_nominal_voltage[3].SetPolar(bus[indexer].volt_base,0.0);
04992 adjust_temp_nominal_voltage[4].SetPolar(bus[indexer].volt_base,-2.0*PI/3.0);
04993 adjust_temp_nominal_voltage[5].SetPolar(bus[indexer].volt_base,2.0*PI/3.0);
04994
04995
04996 adjust_temp_voltage_mag[3] = bus[indexer].V[0].Mag();
04997 adjust_temp_voltage_mag[4] = bus[indexer].V[1].Mag();
04998 adjust_temp_voltage_mag[5] = bus[indexer].V[2].Mag();
04999
05000
05001 if ((bus[indexer].I[0] != 0.0) && (adjust_temp_voltage_mag[3] != 0.0))
05002 {
05003
05004 adjusted_constant_current[0] = ~(adjust_temp_nominal_voltage[3] * ~bus[indexer].I[0] * adjust_temp_voltage_mag[3] / (bus[indexer].V[0] * adjust_nominal_voltage_val));
05005 }
05006 else
05007 {
05008 adjusted_constant_current[0] = complex(0.0,0.0);
05009 }
05010
05011
05012 if ((bus[indexer].I[1] != 0.0) && (adjust_temp_voltage_mag[4] != 0.0))
05013 {
05014
05015 adjusted_constant_current[1] = ~(adjust_temp_nominal_voltage[4] * ~bus[indexer].I[1] * adjust_temp_voltage_mag[4] / (bus[indexer].V[1] * adjust_nominal_voltage_val));
05016 }
05017 else
05018 {
05019 adjusted_constant_current[1] = complex(0.0,0.0);
05020 }
05021
05022
05023 if ((bus[indexer].I[2] != 0.0) && (adjust_temp_voltage_mag[5] != 0.0))
05024 {
05025
05026 adjusted_constant_current[2] = ~(adjust_temp_nominal_voltage[5] * ~bus[indexer].I[2] * adjust_temp_voltage_mag[5] / (bus[indexer].V[2] * adjust_nominal_voltage_val));
05027 }
05028 else
05029 {
05030 adjusted_constant_current[2] = complex(0.0,0.0);
05031 }
05032
05033 if (bus[indexer].prerot_I[0] != 0.0)
05034 adjusted_constant_current[0] += bus[indexer].prerot_I[0];
05035
05036 if (bus[indexer].prerot_I[1] != 0.0)
05037 adjusted_constant_current[1] += bus[indexer].prerot_I[1];
05038
05039 if (bus[indexer].prerot_I[2] != 0.0)
05040 adjusted_constant_current[2] += bus[indexer].prerot_I[2];
05041
05042
05043 if ((bus[indexer].phases & 0x10) == 0x10)
05044 {
05045
05046 adjust_nominal_voltage_val = bus[indexer].volt_base * sqrt(3.0);
05047
05048
05049 adjust_temp_nominal_voltage[0].SetPolar(adjust_nominal_voltage_val,PI/6.0);
05050 adjust_temp_nominal_voltage[1].SetPolar(adjust_nominal_voltage_val,-1.0*PI/2.0);
05051 adjust_temp_nominal_voltage[2].SetPolar(adjust_nominal_voltage_val,5.0*PI/6.0);
05052
05053
05054 voltageDel[0] = bus[indexer].V[0] - bus[indexer].V[1];
05055 voltageDel[1] = bus[indexer].V[1] - bus[indexer].V[2];
05056 voltageDel[2] = bus[indexer].V[2] - bus[indexer].V[0];
05057
05058
05059 adjust_temp_voltage_mag[0] = voltageDel[0].Mag();
05060 adjust_temp_voltage_mag[1] = voltageDel[1].Mag();
05061 adjust_temp_voltage_mag[2] = voltageDel[2].Mag();
05062
05063
05064 if ((bus[indexer].extra_var[6] != 0.0) && (adjust_temp_voltage_mag[0] != 0.0))
05065 {
05066
05067 adjusted_constant_current[3] = ~(adjust_temp_nominal_voltage[0] * ~bus[indexer].extra_var[6] * adjust_temp_voltage_mag[0] / (voltageDel[0] * adjust_nominal_voltage_val));
05068 }
05069 else
05070 {
05071 adjusted_constant_current[3] = complex(0.0,0.0);
05072 }
05073
05074
05075 if ((bus[indexer].extra_var[7] != 0.0) && (adjust_temp_voltage_mag[1] != 0.0))
05076 {
05077
05078 adjusted_constant_current[4] = ~(adjust_temp_nominal_voltage[1] * ~bus[indexer].extra_var[7] * adjust_temp_voltage_mag[1] / (voltageDel[1] * adjust_nominal_voltage_val));
05079 }
05080 else
05081 {
05082 adjusted_constant_current[4] = complex(0.0,0.0);
05083 }
05084
05085
05086 if ((bus[indexer].extra_var[8] != 0.0) && (adjust_temp_voltage_mag[2] != 0.0))
05087 {
05088
05089 adjusted_constant_current[5] = ~(adjust_temp_nominal_voltage[2] * ~bus[indexer].extra_var[8] * adjust_temp_voltage_mag[2] / (voltageDel[2] * adjust_nominal_voltage_val));
05090 }
05091 else
05092 {
05093 adjusted_constant_current[5] = complex(0.0,0.0);
05094 }
05095 }
05096 else
05097 {
05098
05099 adjusted_constant_current[3] = complex(0.0,0.0);
05100 adjusted_constant_current[4] = complex(0.0,0.0);
05101 adjusted_constant_current[5] = complex(0.0,0.0);
05102 }
05103 }
05104 else
05105 {
05106 adjusted_constant_current[0] = bus[indexer].I[0];
05107 adjusted_constant_current[1] = bus[indexer].I[1];
05108 adjusted_constant_current[2] = bus[indexer].I[2];
05109
05110
05111 if ((bus[indexer].phases & 0x10) == 0x10)
05112 {
05113
05114 adjusted_constant_current[3] = bus[indexer].extra_var[6];
05115 adjusted_constant_current[4] = bus[indexer].extra_var[7];
05116 adjusted_constant_current[5] = bus[indexer].extra_var[8];
05117 }
05118 else
05119 {
05120 adjusted_constant_current[3] = complex(0.0,0.0);
05121 adjusted_constant_current[4] = complex(0.0,0.0);
05122 adjusted_constant_current[5] = complex(0.0,0.0);
05123 }
05124 }
05125
05126
05127 temp_index = -1;
05128 temp_index_b = -1;
05129
05130 if ((bus[indexer].phases & 0x10) == 0x10)
05131 {
05132
05133 if ((bus[indexer].phases & 0x06) == 0x06)
05134 {
05135
05136 voltageDel[0] = bus[indexer].V[0] - bus[indexer].V[1];
05137
05138
05139 delta_current[0] = (voltageDel[0] == 0) ? 0 : ~(bus[indexer].extra_var[0]/voltageDel[0]);
05140
05141
05142 delta_current[0] += voltageDel[0] * (bus[indexer].extra_var[3]);
05143 }
05144 else
05145 {
05146
05147 voltageDel[0] = complex(0.0,0.0);
05148 delta_current[0] = complex(0.0,0.0);
05149 }
05150
05151
05152 if ((bus[indexer].phases & 0x03) == 0x03)
05153 {
05154
05155 voltageDel[1] = bus[indexer].V[1] - bus[indexer].V[2];
05156
05157
05158 delta_current[1] = (voltageDel[1] == 0) ? 0 : ~(bus[indexer].extra_var[1]/voltageDel[1]);
05159
05160
05161 delta_current[1] += voltageDel[1] * (bus[indexer].extra_var[4]);
05162 }
05163 else
05164 {
05165
05166 voltageDel[1] = complex(0.0,0.0);
05167 delta_current[1] = complex(0.0,0.0);
05168 }
05169
05170
05171 if ((bus[indexer].phases & 0x05) == 0x05)
05172 {
05173
05174 voltageDel[2] = bus[indexer].V[2] - bus[indexer].V[0];
05175
05176
05177 delta_current[2] = (voltageDel[2] == 0) ? 0 : ~(bus[indexer].extra_var[2]/voltageDel[2]);
05178
05179
05180 delta_current[2] += voltageDel[2] * (bus[indexer].extra_var[5]);
05181 }
05182 else
05183 {
05184
05185 voltageDel[2] = complex(0.0,0.0);
05186 delta_current[2] = complex(0.0,0.0);
05187 }
05188
05189
05190 undeltacurr[0]=(adjusted_constant_current[3]+delta_current[0])-(adjusted_constant_current[5]+delta_current[2]);
05191 undeltacurr[1]=(adjusted_constant_current[4]+delta_current[1])-(adjusted_constant_current[3]+delta_current[0]);
05192 undeltacurr[2]=(adjusted_constant_current[5]+delta_current[2])-(adjusted_constant_current[4]+delta_current[1]);
05193 }
05194 else
05195 {
05196 undeltacurr[0] = undeltacurr[1] = undeltacurr[2] = complex(0.0,0.0);
05197 }
05198
05199 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
05200 {
05201 switch(bus[indexer].phases & 0x07) {
05202 case 0x01:
05203 {
05204 temp_index=0;
05205 temp_index_b=2;
05206 break;
05207 }
05208 case 0x02:
05209 {
05210 temp_index=0;
05211 temp_index_b=1;
05212 break;
05213 }
05214 case 0x03:
05215 {
05216 if (jindex==0)
05217 {
05218 temp_index=0;
05219 temp_index_b=1;
05220 }
05221 else
05222 {
05223 temp_index=1;
05224 temp_index_b=2;
05225 }
05226 break;
05227 }
05228 case 0x04:
05229 {
05230 temp_index=0;
05231 temp_index_b=0;
05232 break;
05233 }
05234 case 0x05:
05235 {
05236 if (jindex==0)
05237 {
05238 temp_index=0;
05239 temp_index_b=0;
05240 }
05241 else
05242 {
05243 temp_index=1;
05244 temp_index_b=2;
05245 }
05246 break;
05247 }
05248 case 0x06:
05249 case 0x07:
05250 {
05251 temp_index=jindex;
05252 temp_index_b=jindex;
05253 break;
05254 }
05255 default:
05256 break;
05257 }
05258
05259 if (jacobian_pass == false)
05260 {
05261 if ((temp_index==-1) || (temp_index_b==-1))
05262 {
05263 GL_THROW("NR: A scheduled power update element failed.");
05264
05265
05266
05267
05268
05269 }
05270
05271
05272 tempPbus = (bus[indexer].S[temp_index_b]).Re();
05273 tempPbus += (adjusted_constant_current[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Re() + (adjusted_constant_current[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Im();
05274 tempPbus += (undeltacurr[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Re() + (undeltacurr[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Im();
05275 tempPbus += (bus[indexer].Y[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Re() + (bus[indexer].Y[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Im();
05276 bus[indexer].PL[temp_index] = tempPbus;
05277
05278
05279 tempQbus = (bus[indexer].S[temp_index_b]).Im();
05280 tempQbus += (adjusted_constant_current[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Im() - (adjusted_constant_current[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Re();
05281 tempQbus += (undeltacurr[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Im() - (undeltacurr[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Re();
05282 tempQbus += -(bus[indexer].Y[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Im() - (bus[indexer].Y[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Re();
05283 bus[indexer].QL[temp_index] = tempQbus;
05284 }
05285 else
05286 {
05287 if ((temp_index==-1) || (temp_index_b==-1))
05288 {
05289 GL_THROW("NR: A Jacobian update element failed.");
05290
05291
05292
05293
05294
05295 }
05296
05297 if ((bus[indexer].V[temp_index_b]).Mag()!=0)
05298 {
05299 bus[indexer].Jacob_A[temp_index] = ((bus[indexer].S[temp_index_b]).Im() * (pow((bus[indexer].V[temp_index_b]).Re(),2) - pow((bus[indexer].V[temp_index_b]).Im(),2)) - 2*(bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(bus[indexer].S[temp_index_b]).Re())/pow((bus[indexer].V[temp_index_b]).Mag(),4);
05300 bus[indexer].Jacob_A[temp_index] += ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(adjusted_constant_current[temp_index_b]).Re() + (adjusted_constant_current[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3) + (bus[indexer].Y[temp_index_b]).Im();
05301 bus[indexer].Jacob_A[temp_index] += ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Re() + (undeltacurr[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05302
05303 bus[indexer].Jacob_B[temp_index] = ((bus[indexer].S[temp_index_b]).Re() * (pow((bus[indexer].V[temp_index_b]).Re(),2) - pow((bus[indexer].V[temp_index_b]).Im(),2)) + 2*(bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(bus[indexer].S[temp_index_b]).Im())/pow((bus[indexer].V[temp_index_b]).Mag(),4);
05304 bus[indexer].Jacob_B[temp_index] += -((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(adjusted_constant_current[temp_index_b]).Im() + (adjusted_constant_current[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3) - (bus[indexer].Y[temp_index_b]).Re();
05305 bus[indexer].Jacob_B[temp_index] += -((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Im() + (undeltacurr[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05306
05307 bus[indexer].Jacob_C[temp_index] = ((bus[indexer].S[temp_index_b]).Re() * (pow((bus[indexer].V[temp_index_b]).Im(),2) - pow((bus[indexer].V[temp_index_b]).Re(),2)) - 2*(bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(bus[indexer].S[temp_index_b]).Im())/pow((bus[indexer].V[temp_index_b]).Mag(),4);
05308 bus[indexer].Jacob_C[temp_index] +=((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(adjusted_constant_current[temp_index_b]).Im() - (adjusted_constant_current[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3) - (bus[indexer].Y[temp_index_b]).Re();
05309 bus[indexer].Jacob_C[temp_index] +=((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Im() - (undeltacurr[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05310
05311 bus[indexer].Jacob_D[temp_index] = ((bus[indexer].S[temp_index_b]).Im() * (pow((bus[indexer].V[temp_index_b]).Re(),2) - pow((bus[indexer].V[temp_index_b]).Im(),2)) - 2*(bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(bus[indexer].S[temp_index_b]).Re())/pow((bus[indexer].V[temp_index_b]).Mag(),4);
05312 bus[indexer].Jacob_D[temp_index] += ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(adjusted_constant_current[temp_index_b]).Re() - (adjusted_constant_current[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3) - (bus[indexer].Y[temp_index_b]).Im();
05313 bus[indexer].Jacob_D[temp_index] += ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Re() - (undeltacurr[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05314
05315 }
05316 else
05317 {
05318 bus[indexer].Jacob_A[temp_index]= (bus[indexer].Y[temp_index_b]).Im() - 1e-4;
05319 bus[indexer].Jacob_B[temp_index]= -(bus[indexer].Y[temp_index_b]).Re() - 1e-4;
05320 bus[indexer].Jacob_C[temp_index]= -(bus[indexer].Y[temp_index_b]).Re() - 1e-4;
05321 bus[indexer].Jacob_D[temp_index]= -(bus[indexer].Y[temp_index_b]).Im() - 1e-4;
05322 }
05323 }
05324 }
05325 }
05326
05327
05328 if ((bus[indexer].phases & 0x80) != 0x80)
05329 {
05330
05331 if ((bus[indexer].phases & 0x06) == 0x06)
05332 {
05333
05334 voltageDel[0] = bus[indexer].V[0] - bus[indexer].V[1];
05335
05336
05337 delta_current[0] = (voltageDel[0] == 0) ? 0 : ~(bus[indexer].S_dy[0]/voltageDel[0]);
05338
05339
05340 delta_current[0] += voltageDel[0] * (bus[indexer].Y_dy[0]);
05341
05342 }
05343 else
05344 {
05345
05346 voltageDel[0] = complex(0.0,0.0);
05347 delta_current[0] = complex(0.0,0.0);
05348 }
05349
05350 if ((bus[indexer].phases & 0x03) == 0x03)
05351 {
05352
05353 voltageDel[1] = bus[indexer].V[1] - bus[indexer].V[2];
05354
05355
05356 delta_current[1] = (voltageDel[1] == 0) ? 0 : ~(bus[indexer].S_dy[1]/voltageDel[1]);
05357
05358
05359 delta_current[1] += voltageDel[1] * (bus[indexer].Y_dy[1]);
05360
05361 }
05362 else
05363 {
05364
05365 voltageDel[1] = complex(0.0,0.0);
05366 delta_current[1] = complex(0.0,0.0);
05367 }
05368
05369 if ((bus[indexer].phases & 0x05) == 0x05)
05370 {
05371
05372 voltageDel[2] = bus[indexer].V[2] - bus[indexer].V[0];
05373
05374
05375 delta_current[2] = (voltageDel[2] == 0) ? 0 : ~(bus[indexer].S_dy[2]/voltageDel[2]);
05376
05377
05378 delta_current[2] += voltageDel[2] * (bus[indexer].Y_dy[2]);
05379
05380 }
05381 else
05382 {
05383
05384 voltageDel[2] = complex(0.0,0.0);
05385 delta_current[2] = complex(0.0,0.0);
05386 }
05387
05388
05389 if (*bus[indexer].dynamics_enabled == true)
05390 {
05391
05392 adjust_nominal_voltage_val = bus[indexer].volt_base;
05393 adjust_nominal_voltaged_val = bus[indexer].volt_base * sqrt(3.0);
05394
05395
05396 adjust_temp_nominal_voltage[0].SetPolar(adjust_nominal_voltaged_val,PI/6.0);
05397 adjust_temp_nominal_voltage[1].SetPolar(adjust_nominal_voltaged_val,-1.0*PI/2.0);
05398 adjust_temp_nominal_voltage[2].SetPolar(adjust_nominal_voltaged_val,5.0*PI/6.0);
05399 adjust_temp_nominal_voltage[3].SetPolar(adjust_nominal_voltage_val,0.0);
05400 adjust_temp_nominal_voltage[4].SetPolar(adjust_nominal_voltage_val,-2.0*PI/3.0);
05401 adjust_temp_nominal_voltage[5].SetPolar(adjust_nominal_voltage_val,2.0*PI/3.0);
05402
05403
05404 voltageDel[0] = bus[indexer].V[0] - bus[indexer].V[1];
05405 voltageDel[1] = bus[indexer].V[1] - bus[indexer].V[2];
05406 voltageDel[2] = bus[indexer].V[2] - bus[indexer].V[0];
05407
05408
05409 adjust_temp_voltage_mag[0] = voltageDel[0].Mag();
05410 adjust_temp_voltage_mag[1] = voltageDel[1].Mag();
05411 adjust_temp_voltage_mag[2] = voltageDel[2].Mag();
05412 adjust_temp_voltage_mag[3] = bus[indexer].V[0].Mag();
05413 adjust_temp_voltage_mag[4] = bus[indexer].V[1].Mag();
05414 adjust_temp_voltage_mag[5] = bus[indexer].V[2].Mag();
05415
05416
05417 if ((bus[indexer].I_dy[3] != 0.0) && (adjust_temp_voltage_mag[3] != 0.0))
05418 {
05419
05420 adjusted_constant_current[3] = ~(adjust_temp_nominal_voltage[3] * ~bus[indexer].I_dy[3] * adjust_temp_voltage_mag[3] / (bus[indexer].V[0] * adjust_nominal_voltage_val));
05421 }
05422 else
05423 {
05424 adjusted_constant_current[3] = complex(0.0,0.0);
05425 }
05426
05427
05428 if ((bus[indexer].I_dy[4] != 0.0) && (adjust_temp_voltage_mag[4] != 0.0))
05429 {
05430
05431 adjusted_constant_current[4] = ~(adjust_temp_nominal_voltage[4] * ~bus[indexer].I_dy[4] * adjust_temp_voltage_mag[4] / (bus[indexer].V[1] * adjust_nominal_voltage_val));
05432 }
05433 else
05434 {
05435 adjusted_constant_current[4] = complex(0.0,0.0);
05436 }
05437
05438
05439 if ((bus[indexer].I_dy[5] != 0.0) && (adjust_temp_voltage_mag[5] != 0.0))
05440 {
05441
05442 adjusted_constant_current[5] = ~(adjust_temp_nominal_voltage[5] * ~bus[indexer].I_dy[5] * adjust_temp_voltage_mag[5] / (bus[indexer].V[2] * adjust_nominal_voltage_val));
05443 }
05444 else
05445 {
05446 adjusted_constant_current[5] = complex(0.0,0.0);
05447 }
05448
05449
05450 if ((bus[indexer].I_dy[0] != 0.0) && (adjust_temp_voltage_mag[0] != 0.0))
05451 {
05452
05453 adjusted_constant_current[0] = ~(adjust_temp_nominal_voltage[0] * ~bus[indexer].I_dy[0] * adjust_temp_voltage_mag[0] / (voltageDel[0] * adjust_nominal_voltaged_val));
05454 }
05455 else
05456 {
05457 adjusted_constant_current[0] = complex(0.0,0.0);
05458 }
05459
05460
05461 if ((bus[indexer].I_dy[1] != 0.0) && (adjust_temp_voltage_mag[1] != 0.0))
05462 {
05463
05464 adjusted_constant_current[1] = ~(adjust_temp_nominal_voltage[1] * ~bus[indexer].I_dy[1] * adjust_temp_voltage_mag[1] / (voltageDel[1] * adjust_nominal_voltaged_val));
05465 }
05466 else
05467 {
05468 adjusted_constant_current[1] = complex(0.0,0.0);
05469 }
05470
05471
05472 if ((bus[indexer].I_dy[2] != 0.0) && (adjust_temp_voltage_mag[2] != 0.0))
05473 {
05474
05475 adjusted_constant_current[2] = ~(adjust_temp_nominal_voltage[2] * ~bus[indexer].I_dy[2] * adjust_temp_voltage_mag[2] / (voltageDel[2] * adjust_nominal_voltaged_val));
05476 }
05477 else
05478 {
05479 adjusted_constant_current[2] = complex(0.0,0.0);
05480 }
05481 }
05482 else
05483 {
05484
05485 adjusted_constant_current[0] = bus[indexer].I_dy[0];
05486 adjusted_constant_current[1] = bus[indexer].I_dy[1];
05487 adjusted_constant_current[2] = bus[indexer].I_dy[2];
05488 adjusted_constant_current[3] = bus[indexer].I_dy[3];
05489 adjusted_constant_current[4] = bus[indexer].I_dy[4];
05490 adjusted_constant_current[5] = bus[indexer].I_dy[5];
05491 }
05492
05493
05494
05495
05496 if ((bus[indexer].phases & 0x04) == 0x04)
05497 {
05498 undeltacurr[0]=(adjusted_constant_current[0]+delta_current[0])-(adjusted_constant_current[2]+delta_current[2]);
05499
05500
05501
05502
05503 undeltacurr[0] += (bus[indexer].V[0] == 0) ? 0 : ~(bus[indexer].S_dy[3]/bus[indexer].V[0]);
05504
05505
05506 undeltacurr[0] += bus[indexer].Y_dy[3]*bus[indexer].V[0];
05507
05508
05509 undeltacurr[0] += adjusted_constant_current[3];
05510 }
05511 else
05512 {
05513
05514 undeltacurr[0] = complex(0.0,0.0);
05515 }
05516
05517 if ((bus[indexer].phases & 0x02) == 0x02)
05518 {
05519 undeltacurr[1]=(adjusted_constant_current[1]+delta_current[1])-(adjusted_constant_current[0]+delta_current[0]);
05520
05521
05522
05523
05524 undeltacurr[1] += (bus[indexer].V[1] == 0) ? 0 : ~(bus[indexer].S_dy[4]/bus[indexer].V[1]);
05525
05526
05527 undeltacurr[1] += bus[indexer].Y_dy[4]*bus[indexer].V[1];
05528
05529
05530 undeltacurr[1] += adjusted_constant_current[4];
05531 }
05532 else
05533 {
05534
05535 undeltacurr[1] = complex(0.0,0.0);
05536 }
05537
05538 if ((bus[indexer].phases & 0x01) == 0x01)
05539 {
05540 undeltacurr[2]=(adjusted_constant_current[2]+delta_current[2])-(adjusted_constant_current[1]+delta_current[1]);
05541
05542
05543
05544
05545 undeltacurr[2] += (bus[indexer].V[2] == 0) ? 0 : ~(bus[indexer].S_dy[5]/bus[indexer].V[2]);
05546
05547
05548 undeltacurr[2] += bus[indexer].Y_dy[5]*bus[indexer].V[2];
05549
05550
05551 undeltacurr[2] += adjusted_constant_current[5];
05552 }
05553 else
05554 {
05555
05556 undeltacurr[2] = complex(0.0,0.0);
05557 }
05558
05559
05560
05561 temp_index = -1;
05562 temp_index_b = -1;
05563
05564 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
05565 {
05566 switch(bus[indexer].phases & 0x07) {
05567 case 0x01:
05568 {
05569 temp_index=0;
05570 temp_index_b=2;
05571 break;
05572 }
05573 case 0x02:
05574 {
05575 temp_index=0;
05576 temp_index_b=1;
05577 break;
05578 }
05579 case 0x03:
05580 {
05581 if (jindex==0)
05582 {
05583 temp_index=0;
05584 temp_index_b=1;
05585 }
05586 else
05587 {
05588 temp_index=1;
05589 temp_index_b=2;
05590 }
05591 break;
05592 }
05593 case 0x04:
05594 {
05595 temp_index=0;
05596 temp_index_b=0;
05597 break;
05598 }
05599 case 0x05:
05600 {
05601 if (jindex==0)
05602 {
05603 temp_index=0;
05604 temp_index_b=0;
05605 }
05606 else
05607 {
05608 temp_index=1;
05609 temp_index_b=2;
05610 }
05611 break;
05612 }
05613 case 0x06:
05614 case 0x07:
05615 {
05616 temp_index=jindex;
05617 temp_index_b=jindex;
05618 break;
05619 }
05620 default:
05621 break;
05622 }
05623
05624 if (jacobian_pass == false)
05625 {
05626 if ((temp_index==-1) || (temp_index_b==-1))
05627 {
05628 GL_THROW("NR: A scheduled power update element failed.");
05629
05630 }
05631
05632
05633 tempPbus = (undeltacurr[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Re() + (undeltacurr[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Im();
05634 bus[indexer].PL[temp_index] += tempPbus;
05635
05636
05637 tempQbus = (undeltacurr[temp_index_b]).Re() * (bus[indexer].V[temp_index_b]).Im() - (undeltacurr[temp_index_b]).Im() * (bus[indexer].V[temp_index_b]).Re();
05638 bus[indexer].QL[temp_index] += tempQbus;
05639 }
05640 else
05641 {
05642 if ((temp_index==-1) || (temp_index_b==-1))
05643 {
05644 GL_THROW("NR: A Jacobian update element failed.");
05645
05646 }
05647
05648 if ((bus[indexer].V[temp_index_b]).Mag()!=0)
05649 {
05650
05651 bus[indexer].Jacob_A[temp_index] += ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Re() + (undeltacurr[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05652 bus[indexer].Jacob_B[temp_index] += -((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Im() + (undeltacurr[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05653 bus[indexer].Jacob_C[temp_index] +=((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Im() - (undeltacurr[temp_index_b]).Re() *pow((bus[indexer].V[temp_index_b]).Im(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05654 bus[indexer].Jacob_D[temp_index] += ((bus[indexer].V[temp_index_b]).Re()*(bus[indexer].V[temp_index_b]).Im()*(undeltacurr[temp_index_b]).Re() - (undeltacurr[temp_index_b]).Im() *pow((bus[indexer].V[temp_index_b]).Re(),2))/pow((bus[indexer].V[temp_index_b]).Mag(),3);
05655 }
05656 else
05657 {
05658 bus[indexer].Jacob_A[temp_index] += -1e-4;
05659 bus[indexer].Jacob_B[temp_index] += -1e-4;
05660 bus[indexer].Jacob_C[temp_index] += -1e-4;
05661 bus[indexer].Jacob_D[temp_index] += -1e-4;
05662 }
05663 }
05664 }
05665 }
05666
05667 if (jacobian_pass == true)
05668 {
05669
05670
05671 if (bus[indexer].full_Y_load != NULL)
05672 {
05673
05674
05675 temp_index = -1;
05676 temp_index_b = -1;
05677
05678 for (jindex=0; jindex<powerflow_values->BA_diag[indexer].size; jindex++)
05679 {
05680 switch(bus[indexer].phases & 0x07) {
05681 case 0x01:
05682 {
05683 temp_index=0;
05684 temp_index_b=2;
05685 break;
05686 }
05687 case 0x02:
05688 {
05689 temp_index=0;
05690 temp_index_b=1;
05691 break;
05692 }
05693 case 0x03:
05694 {
05695 if (jindex==0)
05696 {
05697 temp_index=0;
05698 temp_index_b=1;
05699 }
05700 else
05701 {
05702 temp_index=1;
05703 temp_index_b=2;
05704 }
05705 break;
05706 }
05707 case 0x04:
05708 {
05709 temp_index=0;
05710 temp_index_b=0;
05711 break;
05712 }
05713 case 0x05:
05714 {
05715 if (jindex==0)
05716 {
05717 temp_index=0;
05718 temp_index_b=0;
05719 }
05720 else
05721 {
05722 temp_index=1;
05723 temp_index_b=2;
05724 }
05725 break;
05726 }
05727 case 0x06:
05728 case 0x07:
05729 {
05730 temp_index=jindex;
05731 temp_index_b=jindex;
05732 break;
05733 }
05734 default:
05735 break;
05736 }
05737
05738 if ((temp_index==-1) || (temp_index_b==-1))
05739 {
05740 GL_THROW("NR: A Jacobian update element failed.");
05741
05742 }
05743
05744
05745 bus[indexer].Jacob_A[temp_index] += bus[indexer].full_Y_load[temp_index_b].Im();
05746 bus[indexer].Jacob_B[temp_index] += bus[indexer].full_Y_load[temp_index_b].Re();
05747 bus[indexer].Jacob_C[temp_index] += bus[indexer].full_Y_load[temp_index_b].Re();
05748 bus[indexer].Jacob_D[temp_index] -= bus[indexer].full_Y_load[temp_index_b].Im();
05749 }
05750 }
05751 }
05752 }
05753 }
05754 }
05755
05756
05757 STATUS NR_array_structure_free(NR_SOLVER_STRUCT *struct_of_interest,int number_of_islands)
05758 {
05759 int index_val;
05760 SUPERLU_NR_vars *curr_island_superLU_vars;
05761
05762
05763 curr_island_superLU_vars = NULL;
05764
05765
05766 for (index_val=0; index_val<number_of_islands; index_val++)
05767 {
05768
05769 if (struct_of_interest->island_matrix_values[index_val].deltaI_NR != NULL)
05770 gl_free(struct_of_interest->island_matrix_values[index_val].deltaI_NR);
05771
05772 if (struct_of_interest->island_matrix_values[index_val].Y_offdiag_PQ != NULL)
05773 gl_free(struct_of_interest->island_matrix_values[index_val].Y_offdiag_PQ);
05774
05775 if (struct_of_interest->island_matrix_values[index_val].Y_diag_fixed != NULL)
05776 gl_free(struct_of_interest->island_matrix_values[index_val].Y_diag_fixed);
05777
05778 if (struct_of_interest->island_matrix_values[index_val].Y_diag_update != NULL)
05779 gl_free(struct_of_interest->island_matrix_values[index_val].Y_diag_update);
05780
05781 if (struct_of_interest->island_matrix_values[index_val].Y_Amatrix != NULL)
05782 gl_free(struct_of_interest->island_matrix_values[index_val].Y_Amatrix);
05783
05784
05785 if (struct_of_interest->island_matrix_values[index_val].matrices_LU.a_LU != NULL)
05786 gl_free(struct_of_interest->island_matrix_values[index_val].matrices_LU.a_LU);
05787
05788 if (struct_of_interest->island_matrix_values[index_val].matrices_LU.cols_LU != NULL)
05789 gl_free(struct_of_interest->island_matrix_values[index_val].matrices_LU.cols_LU);
05790
05791 if (struct_of_interest->island_matrix_values[index_val].matrices_LU.rhs_LU != NULL)
05792 gl_free(struct_of_interest->island_matrix_values[index_val].matrices_LU.rhs_LU);
05793
05794 if (struct_of_interest->island_matrix_values[index_val].matrices_LU.rows_LU != NULL)
05795 gl_free(struct_of_interest->island_matrix_values[index_val].matrices_LU.rows_LU);
05796
05797 if (struct_of_interest->island_matrix_values[index_val].LU_solver_vars != NULL)
05798 {
05799
05800 if (matrix_solver_method == MM_SUPERLU)
05801 {
05802
05803 curr_island_superLU_vars = (SUPERLU_NR_vars*)struct_of_interest->island_matrix_values[index_val].LU_solver_vars;
05804
05805
05806 if (curr_island_superLU_vars->A_LU.Store != NULL)
05807 gl_free(curr_island_superLU_vars->A_LU.Store);
05808
05809 if (curr_island_superLU_vars->B_LU.Store != NULL)
05810 gl_free(curr_island_superLU_vars->B_LU.Store);
05811
05812 if (curr_island_superLU_vars->perm_c != NULL)
05813 gl_free(curr_island_superLU_vars->perm_c);
05814
05815 if (curr_island_superLU_vars->perm_r != NULL)
05816 gl_free(curr_island_superLU_vars->perm_r);
05817
05818
05819 curr_island_superLU_vars = NULL;
05820
05821
05822 gl_free(struct_of_interest->island_matrix_values[index_val].LU_solver_vars);
05823 }
05824 else if (matrix_solver_method == MM_EXTERN)
05825 {
05826
05827 ((void (*)(void *, bool))(LUSolverFcns.ext_destroy))(struct_of_interest->island_matrix_values[index_val].LU_solver_vars,false);
05828 }
05829
05830 }
05831 }
05832
05833
05834 gl_free(struct_of_interest->island_matrix_values);
05835
05836
05837 struct_of_interest->island_matrix_values = NULL;
05838
05839
05840 gl_free(struct_of_interest->BA_diag);
05841
05842
05843 struct_of_interest->BA_diag = NULL;
05844
05845
05846 return SUCCESS;
05847 }
05848
05849
05850 STATUS NR_array_structure_allocate(NR_SOLVER_STRUCT *struct_of_interest,int number_of_islands)
05851 {
05852 int index_val;
05853
05854
05855 if (struct_of_interest->island_matrix_values != NULL)
05856 {
05857 gl_error("solver_NR: Attempted to allocate over an existing NR solver variable array!");
05858
05859
05860
05861
05862
05863
05864 return FAILED;
05865 }
05866
05867
05868 struct_of_interest->island_matrix_values = (NR_MATRIX_CONSTRUCTION *)gl_malloc(number_of_islands*sizeof(NR_MATRIX_CONSTRUCTION));
05869
05870
05871 if (struct_of_interest->island_matrix_values == NULL)
05872 {
05873 gl_error("solver_NR: An attempt to allocate a variable array for the NR solver failed!");
05874
05875
05876
05877
05878
05879 return FAILED;
05880 }
05881
05882
05883 for (index_val=0; index_val<number_of_islands; index_val++)
05884 {
05885
05886 struct_of_interest->island_matrix_values[index_val].deltaI_NR = NULL;
05887 struct_of_interest->island_matrix_values[index_val].size_offdiag_PQ = 0;
05888 struct_of_interest->island_matrix_values[index_val].size_diag_fixed = 0;
05889 struct_of_interest->island_matrix_values[index_val].total_variables = 0;
05890 struct_of_interest->island_matrix_values[index_val].size_diag_update = 0;
05891 struct_of_interest->island_matrix_values[index_val].size_Amatrix = 0;
05892 struct_of_interest->island_matrix_values[index_val].max_size_offdiag_PQ = 0;
05893 struct_of_interest->island_matrix_values[index_val].max_size_diag_fixed = 0;
05894 struct_of_interest->island_matrix_values[index_val].max_total_variables = 0;
05895 struct_of_interest->island_matrix_values[index_val].max_size_diag_update = 0;
05896 struct_of_interest->island_matrix_values[index_val].prev_m = 0;
05897 struct_of_interest->island_matrix_values[index_val].index_count = 0;
05898 struct_of_interest->island_matrix_values[index_val].bus_count = 0;
05899 struct_of_interest->island_matrix_values[index_val].indexer = 0;
05900 struct_of_interest->island_matrix_values[index_val].NR_realloc_needed = false;
05901 struct_of_interest->island_matrix_values[index_val].Y_offdiag_PQ = NULL;
05902 struct_of_interest->island_matrix_values[index_val].Y_diag_fixed = NULL;
05903 struct_of_interest->island_matrix_values[index_val].Y_diag_update = NULL;
05904 struct_of_interest->island_matrix_values[index_val].Y_Amatrix = NULL;
05905
05906
05907 struct_of_interest->island_matrix_values[index_val].matrices_LU.a_LU = NULL;
05908 struct_of_interest->island_matrix_values[index_val].matrices_LU.cols_LU = NULL;
05909 struct_of_interest->island_matrix_values[index_val].matrices_LU.rhs_LU = NULL;
05910 struct_of_interest->island_matrix_values[index_val].matrices_LU.rows_LU = NULL;
05911
05912
05913 struct_of_interest->island_matrix_values[index_val].LU_solver_vars = NULL;
05914 struct_of_interest->island_matrix_values[index_val].iteration_count = 0;
05915 struct_of_interest->island_matrix_values[index_val].new_iteration_required = false;
05916 struct_of_interest->island_matrix_values[index_val].swing_converged = false;
05917 struct_of_interest->island_matrix_values[index_val].swing_is_a_swing = false;
05918 struct_of_interest->island_matrix_values[index_val].SaturationMismatchPresent = false;
05919 struct_of_interest->island_matrix_values[index_val].solver_info = -1;
05920 struct_of_interest->island_matrix_values[index_val].return_code = -1;
05921 struct_of_interest->island_matrix_values[index_val].max_mismatch_converge = 0.0;
05922 }
05923
05924
05925 struct_of_interest->BA_diag = NULL;
05926
05927
05928 return SUCCESS;
05929 }