/***************************************************************************************** * ENS 1.0 - Electric Network Simulator * * Author: * * Lluís Ros (c) 1995-2007 * Institut de Robōtica i Informātica Industrial (CSIC-UPC) * Llorens Artigas 4-6, 2a planta * 08028 Barcelona * Web: http://www-iri.upc.es/people/ros * E-mail: llros@iri.upc.es * * What is ENS? * * ENS is a library of functions to simulate the state of an electric * distribution network. The simulator implements, with some variations, * the sweeping/compensation method by Luo and Semlyen, published in * * G. X. Luo and A. Semlyen. "Efficient load flow for large weakly * meshed networks", IEEE Trans. on Power Systems, Vol. 5, No. 4, * November 1990. * * The network state is computed from the following input: * * - The complex voltage at source buses. * - The complex consumed or injected power at all PQ buses * (the sum of P and Q for all bus loads, generations, or * shunt capacitors connected to the bus). * - The voltage magnitude and real power at PV buses. * - The complex impedance of all branches. * * All complex bus voltages and complex branch currents are computed * as output. Cycles (i.e., loops, or meshes) are allowed to exist * in the network, in any amount. Null-impedance cycles may also be * present. Shunt capacitors with known PQ values may easily be modelled as * * Associated files: * * ens.c, ens.h - functions to simulate the network * glr.c, glr.h - functions to handle the network graph * llr.c, llr.h - functions to handle lists * * Conditions of use: * * The following software may be used by anyone (hereafter the "user") agreeing with the * following conditions: * * 1. The user will only use the software for nonprofit, nonmilitary purposes. * Any commercial use of this software requires the obtention of an appropriate * software license jointly issued by the author and the Institut de Robōtica * i Informātica Industrial (UPC/CSIC). Should you be interested in such license, * please contact the author at the address above. * 2. On any publication reporting results using this software, the user commits * to make a reference to the software, as follows: * * L. Ros, "ENS: An Electric Network Simulator for Distribution * Networks". Institut de Robōtica i Informātica Industrial, Barcelona, * Spain. Available through http://www-iri.upc.es/people/ros * * 3. The user will report the software's author about any trouble he meets in the use of * this software. The author will remove eventual bugs at his earliest convenience. * 4. The author assumes no responsibility for any trouble or damage caused by the * user when using this software for any purpose. * 5. The user will notify the author (to the address given above) that she will use the * software. In her notification the user is asked to please identify herself, to state * she accepts these conditions, and briefly report the purpose for which the software * will be used. * * The functions are organized into the following sections: * * Section 1 - Network constructors * Functions to read an electric network from file * Section 2 - Network destructors * Functions to destroy an electric network * Section 3 - Network simulations * Functions to compute a power flow (for the moment) * Section 4 - Network manoeuvers * Functions to open/close switches and locally recompute the network state * Section 5 - Miscellaneous functions * Miscellaneous functions * *********************************************************************************************/ #include "ens.h" /* * Macros to define the verbosity level: * #define TURN_ON_PRINT1 activates the printing of "essential" output * #define TURN_ON_PRINT2 activates the printing of more detailed output */ #define TURN_ON_PRINT1 #include "debug_macros.h" /* * Private variables */ static double max_discrepancy; /* * Private function prototypes */ static boolean is_one_of(unsigned int num, unsigned int* numbers, int nnumbers); /****************************************************************************************************** * * Section 1 - Network constructors * ******************************************************************************************************/ status init_network(graph* p_N) /* * Initializes a network */ { return(init_graph(p_N)); } /*----------------------------------------------------------------------------------------------------*/ void load_network(FILE* fd, graph* p_N, boolean allareas, area_id_type* areas, int nareas) /* * Loads a network from a file in ENS format. Assumes that all bus fields in * the file appear before the branch fields. Only the network areas specified in * areas[0,...,nareas-1] are loaded. If allareas==TRUE, then all areas are loaded. */ { char file_line[200]; char element[20]; /* Explore all lines of the network file. If areas!=NULL, only load buses specified in "areas" */ while(fgets(file_line, 200, fd)!=NULL) { /* Avoid comments and empty lines */ if(file_line[0]!='#' && file_line[0]!='\n') { /* Parse first element of the line */ sscanf(file_line, "%s", element); /* Check element type */ if(strcmp(element,"BUS")==0) { /* Element is a bus */ load_bus(file_line, p_N, allareas, areas, nareas); } else if(strcmp(element,"BRANCH")==0) { /* Element is a branch */ load_branch(file_line, p_N, allareas, areas, nareas); } else { /* Unknown element */ ens_error("Error in load_network():: invalid line in file\n"); } } } } /*----------------------------------------------------------------------------------------------------*/ void load_bus(char* file_line, graph* p_N, boolean allareas, area_id_type* areas, int nareas) /* * Loads, from a file, a bus into the graph */ { /* Temporal vars to store a parsed bus record */ char vertex_bus_tag[10]; vertex_id_type vertex_id; area_id_type vertex_area_id; int vertex_type; double vertex_inj_pwr_r; double vertex_inj_pwr_i; double vertex_voltage_r; double vertex_voltage_i; double vertex_base_power; double vertex_base_voltage; /* Other vars */ specific_vertex_data* p_svd; /* Parse the bus record */ sscanf(file_line, "%s %u %u %d %lf %lf %lf %lf %lf %lf", vertex_bus_tag, &vertex_id, &vertex_area_id, &vertex_type, &vertex_inj_pwr_r, &vertex_inj_pwr_i, &vertex_voltage_r, &vertex_voltage_i, &vertex_base_power, &vertex_base_voltage ); /* If bus has to be loaded... */ if(allareas==TRUE || is_one_of(vertex_area_id,areas,nareas)) { /* Allocate specific vertex data struct and fill it in */ p_svd=(specific_vertex_data*)malloc(sizeof(specific_vertex_data)); if(p_svd==NULL) ens_error("Error:: load_bus() couldn't allocate memory for p_svd"); p_svd->vertex_id = vertex_id; p_svd->vertex_area = vertex_area_id; p_svd->vertex_type = vertex_type; p_svd->vertex_state = UNKNOWN; p_svd->vertex_injected_pwr.r = vertex_inj_pwr_r; p_svd->vertex_injected_pwr.i = vertex_inj_pwr_i; p_svd->vertex_accumulated_pwr = Complex(0.,0.); p_svd->vertex_voltage.r = vertex_voltage_r; p_svd->vertex_voltage.i = vertex_voltage_i; p_svd->vertex_base_power = vertex_base_power; p_svd->vertex_base_voltage = vertex_base_voltage; /* Add vertex to graph */ if(add_vertex(p_N, vertex_id, (generic_ptr)p_svd)==ERROR) { ens_error("Error:: load_network_areas() couldn't add vertex"); } else { PRINT2(("Loaded bus: %u %u %d %.4e %.4e %.4e %.4e %.4e %.4e\n", vertex_id, vertex_area_id, vertex_type, vertex_inj_pwr_r, vertex_inj_pwr_i, vertex_voltage_r, vertex_voltage_i, vertex_base_power, vertex_base_voltage)); } } } /*----------------------------------------------------------------------------------------------------*/ void load_branch(char* file_line, graph* p_N, boolean allareas, area_id_type* areas, int nareas) /* * Loads, from a file, a network branch into the bus */ { /* Temporal vars to store a parsed branch record */ char edge_branch_tag[10]; edge_id_type edge_id; int edge_type; vertex_id_type edge_vertex1; vertex_id_type edge_vertex2; area_id_type edge_area1; area_id_type edge_area2; int edge_state; int edge_state_nominal; double edge_impedance_r; double edge_impedance_i; double edge_intensity_r; double edge_intensity_i; double edge_max_intensity; /* Other vars */ specific_edge_data* p_sed; status exit_code; /* Parse the bus record */ sscanf(file_line, "%s %u %d %u %u %u %u %d %d %lf %lf %lf %lf %lf", edge_branch_tag, &edge_id, &edge_type, &edge_vertex1, &edge_vertex2, &edge_area1, &edge_area2, &edge_state, &edge_state_nominal, &edge_impedance_r, &edge_impedance_i, &edge_intensity_r, &edge_intensity_i, &edge_max_intensity ); /* If the branch has to be loaded... */ if(allareas==TRUE || (is_one_of(edge_area1,areas,nareas) && is_one_of(edge_area2,areas,nareas))) { /* Allocate specific edge data struct */ p_sed=(specific_edge_data*)malloc(sizeof(specific_edge_data)); if(p_sed==NULL) ens_error("Error:: load_branch() couldn't allocate memory for p_sed"); /* Fill it in */ p_sed->edge_id = edge_id; p_sed->edge_type = edge_type; p_sed->edge_vertex1 = edge_vertex1; p_sed->edge_vertex2 = edge_vertex2; p_sed->edge_state = edge_state; p_sed->edge_state_nominal = edge_state_nominal; p_sed->edge_impedance.r = edge_impedance_r; p_sed->edge_impedance.i = edge_impedance_i; p_sed->edge_intensity.r = edge_intensity_r; p_sed->edge_intensity.i = edge_intensity_i; p_sed->edge_max_intensity = edge_max_intensity; /* Add edge to graph */ exit_code = add_edge(*p_N, edge_vertex1, edge_vertex2, (generic_ptr)p_sed); if(exit_code == ERROR) { ens_error("Error: load_branch() couldn't add_edge to graph"); } else { PRINT2(("Loaded branch: %u %d %u %u %u %u %d %d %.4e %.4e %.4e %.4e %.4e\n", edge_id, edge_type, edge_vertex1, edge_vertex2, edge_area1, edge_area2, edge_state, edge_state_nominal, edge_impedance_r, edge_impedance_i, edge_intensity_r, edge_intensity_i, edge_max_intensity )); } } } /***************************************************************************************************** * * Section 3 - Network simulation * *****************************************************************************************************/ void power_flow_network(graph N) /* * Computes the electric state of the whole network */ { list vertices; graph_vertex vtx; /* Set up a list with all buses of the network */ init_list(&vertices); for(vtx=N;vtx!=NULL;vtx=NEXT(vtx)) insert(&vertices,(generic_ptr)vtx); /* Compute the power flow */ power_flow_subnetwork(vertices); /* Free list of buses */ destroy_list(&vertices,NULL); } /*----------------------------------------------------------------------------------------------------*/ void power_flow_subnetwork(list vertices) /* * Computes the electric state of the subnetwork reachable from a set of vertices. * Warning: the function uses max_discrepancy, which is a global variable. */ { int iter; clock_t tbef, taft; status err; list feeders; list chords; int dim; double** Z; int* indx; double d; boolean network_has_cycles; /* Initializations */ err=init_list(&feeders); if(err==ERROR) ens_error("Error::failed to initialize list of feeders in power_flow()"); err=init_list(&chords); if(err==ERROR) ens_error("Error::failed to initialize list of chords in power_flow()"); /* Check input */ if(vertices==NULL) ens_error("Error:: Empty network received in power_flow()"); /* Reset subnetwork (marks reachable subnetwork as "de-energized") */ reset_subnetwork(vertices); /* Get the list of feeders connected to vertices */ find_feeders(vertices, &feeders, branch_state_is_closed); /* Determine the chords (breakpoints) and set parent pointers for all vertices */ find_chords(feeders,&chords,branch_state_is_closed); /* Mark the positive and negative paths determining the cycles */ mark_negative_positive_paths(chords); /* If applicable, compute the sensitivity matrix, LU decompose it, and open the cycles */ dim = count_non_null_impedance_cycles(chords); network_has_cycles = (dim > 0); if(network_has_cycles) { /* Open the chord (breakpoint) edges */ open_chords(chords); /* Compute sensitivity matrix */ Z = dmatrix(1,2*dim,1,2*dim); compute_sensitivity_matrix(chords,Z,dim); PRINT1(("Sensitivity matrix is:\n")); print_dmatrix(Z,1,2*dim,1,2*dim); /* LU decompose it */ indx = ivector(1,2*dim); dludcmp(Z,2*dim,indx,&d); PRINT1(("LU decomposition is:\n")); print_dmatrix(Z,1,2*dim,1,2*dim); } /* Init. voltages and intensities to 1+j0 pu and 0+j0 pu */ init_subnetwork(feeders); /* Sweeps */ tbef = clock(); iter = 0; do { iter ++; max_discrepancy = 0.; sweep_back_sv(feeders); sweep_forw_sv(feeders); if(network_has_cycles) { inject_power_corrections(chords,Z,dim,indx); } PRINT1(("Iteration %d: max_discrepancy = %.2e\n",iter,max_discrepancy)); } while (iter < MAX_ITER && max_discrepancy > MAX_ERROR); taft = clock(); /* Diagnose */ if(iter==MAX_ITER) PRINT1(("Error:: power_flow exceeded the maximum number of iterations\n")); else PRINT1(("Power flow completed:\n")); /* Iterations performed and CPU time used */ PRINT1(("\tIterations done: %d\n",iter)); PRINT1(("\tMax. discrepancy: %.2e\n", max_discrepancy)); PRINT1(("\tEllapsed time: %.3f sec\n", (double)((taft-tbef)/CLOCKS_PER_SEC))); /* Reset max_discrepancy for later use */ max_discrepancy = 0.; if(network_has_cycles) { /* Restore the chord states */ close_chords(chords); /* Unlabel positive/negative paths */ unmark_negative_positive_paths(chords); /* Free dynamic matrices and vectors */ free_dmatrix(Z,1,2*dim,1,2*dim); free_ivector(indx,1,2*dim); } /* Free dynamic memory used */ destroy_list(&feeders,NULL); destroy_list(&chords,free); } /*----------------------------------------------------------------------------------------------------*/ int count_non_null_impedance_cycles(list chords) { list ch; int count; dcomplex z; count = 0; for(ch=chords;ch!=NULL;ch=NEXT(ch)) { z = self_impedance(ch); if(non_tiny_complex(z)) count++; } return(count); } /*----------------------------------------------------------------------------------------------------*/ void inject_power_corrections(list chords, double** Z, int dim, int* indx) /* * Given the breakpoints in "chords", and their LU decomposed impedance matrix * Z[1..dim][1..dim], with associated index vector indx[1..dim], this function * computes the power injections required at the breakpoints, according to the * compensation method by Luo and Semlyen. These are introduced in the graph * for later use by the backward/forward sweepings. */ { list ch; /* Current chord */ graph_vertex kp; /* Pointer to vertex k' of the chord */ graph_vertex k; /* Pointer to vertex k of the chord */ double* dv; /* Vector of voltage mismatches vk'-vk */ dcomplex v_kp; /* Voltage of vertex k' */ dcomplex v_k; /* Voltage of (virtual) vertex k */ dcomplex v_mis; /* vkp - vk */ dcomplex ds; dcomplex di; int i; /* Voltage mismatches */ dv = dvector(1,2*dim); i=1; for(ch=chords;ch!=NULL;ch=NEXT(ch)) { kp = CHORD_VERTEX1(ch); k = CHORD_VERTEX2(ch); v_kp = VERTEX_VOLTAGE(kp); v_k = VERTEX_VOLTAGE(k); v_mis = Csub(v_kp,v_k); dv[i] = v_mis.r; dv[i+dim] = v_mis.i; i++; } /* Cycle currents (returned in dv itself) */ dlubksb(Z,dim,indx,dv); /* Set chord edge currents to the cycle currents */ i=1; for(ch=chords;ch!=NULL;ch=NEXT(ch)) { EDGE_INTENSITY(CHORD_EDGE(ch)) = Complex(dv[i],dv[i+dim]); PRINT2(("Breakpoint current at edge %u = (%.3e,%.3e)\n", EDGE_ID(CHORD_EDGE(ch)), dv[i], dv[i+dim])); i++; } /* Power corrections */ i=1; for(ch=chords;ch!=NULL;ch=NEXT(ch)) { kp = CHORD_VERTEX1(ch); k = CHORD_VERTEX2(ch); v_kp = VERTEX_VOLTAGE(kp); v_k = VERTEX_VOLTAGE(k); di = Complex(dv[i],dv[i+dim]); ds = Cmul(RCmul(0.5,Cadd(v_k,v_kp)),Conjg(di)); PRINT2(("Breakpoint power correction at edge %u = (%.3e,%.3e)\n", EDGE_ID(CHORD_EDGE(ch)), ds.r, ds.i)); VERTEX_ACUM_PWR(kp) = Csub(VERTEX_ACUM_PWR(kp),ds); VERTEX_ACUM_PWR(k) = Cadd(VERTEX_ACUM_PWR(k), ds); i++; } } /*----------------------------------------------------------------------------------------------------*/ void open_chords(list chords) /* * Marks the edges indicated in chords as "OPEN" */ { list ch; for(ch=chords;ch!=NULL;ch=NEXT(ch)) { EDGE_STATE(CHORD_EDGE(ch)) = OPEN; EDGE_INTENSITY(CHORD_EDGE(ch)) = Complex(0.,0.); } } /*----------------------------------------------------------------------------------------------------*/ void close_chords(list chords) /* * Marks the edges indicated in chords as "CLOSED" */ { list ch; for(ch=chords;ch!=NULL;ch=NEXT(ch)) { EDGE_STATE(CHORD_EDGE(ch)) = CLOSED; } } /*----------------------------------------------------------------------------------------------------*/ void compute_sensitivity_matrix(list chords, double** Z, int dim) /* * Computes the sensitivity matrix associated with a list of chords (breakpoints) * Prerequisites: the associated cycles must have their negative and positive * paths properly labelled. The associated vertices must have proper parent * pointers. */ { list chi, chj; int i, j, labeli, labelj; dcomplex zii, zjj, zij; graph_edge pe; i=1; for(chi=chords; chi!=NULL; chi=NEXT(chi)) { labeli = *((int*)DATA(EDGE_LABELS(CHORD_EDGE(chi)))); zii = self_impedance(chi); /* If the cycle has null impedance, we skip it */ if(non_tiny_complex(zii)) { j=i; for(chj=chi; chj!=NULL; chj=NEXT(chj)) { labelj = *((int*)DATA(EDGE_LABELS(CHORD_EDGE(chj)))); zjj = self_impedance(chj); /* If the cycle has null impedance, we skip it */ if(non_tiny_complex(zjj)) { if(i==j) { Z[i][i] = zii.r; Z[i+dim][i+dim] = zii.r; Z[i][i+dim] = -zii.i; Z[i+dim][i] = zii.i; } else { zij = mutual_impedance(chi,chj); Z[i][j] = zij.r; Z[j][i] = zij.r; Z[i+dim][j+dim] = zij.r; Z[j+dim][i+dim] = zij.r; Z[i][j+dim] = -zij.i; Z[i+dim][j] = zij.i; Z[j][i+dim] = -zij.i; Z[j+dim][i] = zij.i; } j++; } } i++; } } } /*----------------------------------------------------------------------------------------------------*/ dcomplex self_impedance(list ch) /* * Computes the self impedance of the cycle defined by chord "ch" */ { graph_vertex pv; graph_edge pe; dcomplex z; int label; /* Get label of cycle */ label = *((int*)DATA(EDGE_LABELS(CHORD_EDGE(ch)))); /* Count impedance of the chord */ z = EDGE_IMPEDANCE(CHORD_EDGE(ch)); /* Add up impedances of the negative path */ for(pv=CHORD_VERTEX1(ch); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); if(edge_is_labelled(pe,-label)) z = Cadd(z,EDGE_IMPEDANCE(pe)); } /* Add up imepdances of the positive path */ for(pv=CHORD_VERTEX2(ch); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); if(edge_is_labelled(pe,label)) z = Cadd(z,EDGE_IMPEDANCE(pe)); } /* Return the self impedance */ return(z); } /*----------------------------------------------------------------------------------------------------*/ dcomplex mutual_impedance(list chi, list chj) /* * Computes the mutual impedance of the cycles labelled as labeli and labelj */ { graph_vertex pv; graph_edge pe; dcomplex z; int labeli, labelj; /* Get labels of cycles */ labeli = *((int*)DATA(EDGE_LABELS(CHORD_EDGE(chi)))); labelj = *((int*)DATA(EDGE_LABELS(CHORD_EDGE(chj)))); /* Count impedance of the chord */ z = Complex(0.,0.); /* Add common impedances along the negative path of chord chi */ for(pv=CHORD_VERTEX1(chi); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); if(edge_is_labelled(pe,-labeli)) { if(edge_is_labelled(pe,-labelj)) { z = Cadd(z,EDGE_IMPEDANCE(pe)); } else if(edge_is_labelled(pe,labelj)) { z = Csub(z,EDGE_IMPEDANCE(pe)); } } } /* Add common impedances along the positive path of chord chi */ for(pv=CHORD_VERTEX2(chi); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); if(edge_is_labelled(pe,labeli)) { if(edge_is_labelled(pe,labelj)) { z = Cadd(z,EDGE_IMPEDANCE(pe)); } else if(edge_is_labelled(pe,-labelj)) { z = Csub(z,EDGE_IMPEDANCE(pe)); } } } /* Return the mutual impedance */ return(z); } /*----------------------------------------------------------------------------------------------------*/ void mark_negative_positive_paths(list chords) /* * Marks the negative and positive paths associated with a list of chords. */ { graph_vertex pv; graph_edge pe; int label; list ch; label=1; for(ch=chords; ch!=NULL; ch=NEXT(ch)) { /* Label chord */ label_edge(CHORD_EDGE(ch),label); /* Label negative path */ for(pv=CHORD_VERTEX1(ch); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); label_edge(pe,-label); } /* Label positive path */ for(pv=CHORD_VERTEX2(ch); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); label_edge(pe,label); } label++; } } /*----------------------------------------------------------------------------------------------------*/ void unmark_negative_positive_paths(list chords) /* * Unmarks the negative and positive paths associated with a list of chords. */ { graph_vertex pv; graph_edge pe; list ch; for(ch=chords; ch!=NULL; ch=NEXT(ch)) { /* Unlabel chord */ destroy_list(&EDGE_LABELS(CHORD_EDGE(ch)),free); /* Unlabel negative path */ for(pv=CHORD_VERTEX1(ch); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); destroy_list(&EDGE_LABELS(pe),free); } /* Unlabel positive path */ for(pv=CHORD_VERTEX2(ch); PARENT(pv)!=NULL; pv=PARENT(pv)) { pe = find_destin_with_pointer(EMANATING(pv),PARENT(pv)); destroy_list(&EDGE_LABELS(pe),free); } } } /*----------------------------------------------------------------------------------------------------*/ void reset_subnetwork(list vertices) /* * It sets, for all buses in the subnetwork reachable from vertices, their state to DISCONNECTED. * Also, for all reachable buses and branches: * - The bus states to DISCONNECTED * - The bus voltages are set to 0+j0 pu * - The branch intensities are set to 0+j0 pu */ { traverse_graph(BF,vertices,branch_state_is_closed,reset_bus_branch,NULL,NULL,NULL); } /*----------------------------------------------------------------------------------------------------*/ void init_subnetwork(list feeders) /* * It sets, for all buses and branches in the subnetwork reachable from feeders: * - The bus states to CONNECTED * - The bus voltages to 1+j0 pu * - The branch intensities to 0+j0 pu */ { traverse_graph(BF,feeders,branch_state_is_closed,init_bus_branch,NULL,NULL,NULL); } /*----------------------------------------------------------------------------------------------------*/ void sweep_back_sv(list feeders) /* * Performs a backward sweep of powers */ { traverse_graph(BF, feeders, branch_state_is_closed, NULL, calc_acum_power, NULL, NULL); } /*----------------------------------------------------------------------------------------------------*/ void sweep_forw_sv(list feeders) /* * Performs a forward sweep of voltages and intensities */ { traverse_graph(BF, feeders, branch_state_is_closed, calc_voltages_intensities, NULL, NULL, NULL); } /*----------------------------------------------------------------------------------------------------*/ status calc_acum_power(graph_vertex vertex_ptr, graph_vertex parent_vertex_ptr, generic_ptr* input, generic_ptr* output) /* * Computes the accumulated power of vertex_ptr. This is the sum of the * accumulated power of its downstream vertices, plus the power consumed * at incident downstream branches, minus the power injected at vertex_ptr. */ { graph_edge e; /* pointer to traverse the adjacent branches of vertex_ptr */ graph_vertex v; /* pointer to traverse the adjacent vertices of vertex_ptr */ dcomplex Sacum; /* acumulated complex power of vertex_ptr */ dcomplex Sdesc; /* acumulated complex power of a descendant */ dcomplex Se; /* complex power consumed at a descendant branch e */ Sacum = Complex(0., 0.); PRINT2(("Computing accumulated power at vertex %u\n",VERTEX_ID(vertex_ptr))); /* For all adjacent vertices v of vertex_ptr... */ for(e=EMANATING(vertex_ptr);e!=NULL;e=NEXT(e)) { v = VERTEX_PTR(e); /* If the branch e=(vertex_ptr,v) is downstream and active... */ if(v != parent_vertex_ptr && EDGE_STATE(e)==CLOSED) { /* Accumulate the power of "v" with the power consumed at "e" */ Sdesc = VERTEX_ACUM_PWR(v); Se = RCmul( Cabs_sq(VERTEX_ACUM_PWR(v))/Cabs_sq(VERTEX_VOLTAGE(v)), EDGE_IMPEDANCE(e) ); Sacum = Cadd(Sacum,Cadd(Sdesc,Se)); PRINT2(("\tPower of descendent vertex %u = (%lf,%lf)\n", VERTEX_ID(v), Sdesc.r, Sdesc.i)); PRINT2(("\tPower of descendent edge = (%lf,%lf)\n", Se.r, Se.i)); } } /* Subtract the injected power at the node (this also works at a consumer) */ PRINT2(("\tPower injected at vertex %u = (%lf,%lf)\n", VERTEX_ID(vertex_ptr), VERTEX_INJ_PWR(vertex_ptr).r, VERTEX_INJ_PWR(vertex_ptr).i)); VERTEX_ACUM_PWR(vertex_ptr) = Csub(Sacum,VERTEX_INJ_PWR(vertex_ptr)); PRINT2(("\tTotal accumulated power at vertex %u = (%lf,%lf)\n", VERTEX_ID(vertex_ptr), VERTEX_ACUM_PWR(vertex_ptr).r, VERTEX_ACUM_PWR(vertex_ptr).i)); return(OK); } /*----------------------------------------------------------------------------------------------------*/ status calc_voltages_intensities(graph_vertex vertex_ptr, graph_vertex parent_vertex_ptr, generic_ptr* input, generic_ptr* output) /* * Computes the voltage of vertex "vertex_ptr" and the intensity of its preceding branch. * Warning: this function acts on max_discrepancy, which is a global variable */ { graph_edge edge_ptr; dcomplex prev_V; double curr_discrepancy; /* If the parent vertex does not exist, the vertex is a feeder and we do nothing */ if(parent_vertex_ptr!=NULL) { /* Get pointer to edge between the two vertices */ for(edge_ptr=EMANATING(vertex_ptr); edge_ptr!=NULL && VERTEX_PTR(edge_ptr)!=parent_vertex_ptr; edge_ptr=NEXT(edge_ptr)); if(edge_ptr==NULL) ens_error("Error:: couldn't find pointer to edge in calc_voltages_intensities\n"); /* Intensity, estimated as I = (S/Vup)* */ EDGE_INTENSITY(edge_ptr) = Conjg(Cdiv(VERTEX_ACUM_PWR(vertex_ptr), VERTEX_VOLTAGE(parent_vertex_ptr))); PRINT2(("Edge (%u,%u) has intensity (%lf,%lf)\n", VERTEX_ID(vertex_ptr), VERTEX_ID(parent_vertex_ptr), EDGE_INTENSITY(edge_ptr).r, EDGE_INTENSITY(edge_ptr).i)); /* Back-up the previous voltage to compute discrepancies later */ prev_V = VERTEX_VOLTAGE(vertex_ptr); /* Vdown = Vup - Z * (S/Vup)* */ VERTEX_VOLTAGE(vertex_ptr) = Csub( VERTEX_VOLTAGE(parent_vertex_ptr), Cmul( EDGE_IMPEDANCE(edge_ptr), EDGE_INTENSITY(edge_ptr)) ); PRINT2(("Vertex %u has voltage (%lf,%lf)\n", VERTEX_NUMBER(vertex_ptr), VERTEX_VOLTAGE(vertex_ptr).r, VERTEX_VOLTAGE(vertex_ptr).i)); /* Update max_discrepancy; */ curr_discrepancy = Cabs(Csub(VERTEX_VOLTAGE(vertex_ptr),prev_V)); PRINT2(("curr_discrepancy = %lf, max_discrepancy = %lf\n", curr_discrepancy, max_discrepancy)); if(curr_discrepancy > max_discrepancy) max_discrepancy = curr_discrepancy; } return(OK); } /*----------------------------------------------------------------------------------------------------*/ status reset_bus_branch(graph_vertex vertex_ptr, graph_vertex parent_vertex_ptr, generic_ptr* input, generic_ptr* output) /* * It resets the bus and preceding branch. * - The bus state is set to DISCONNECTED * - The bus voltage is set to 0+j0 pu * - The branch intensity is set to 0+j0 pu */ { graph_edge edge_ptr; VERTEX_STATE(vertex_ptr) = DISCONNECTED; VERTEX_VOLTAGE(vertex_ptr) = Complex(0.,0.); /* If the parent vertex does not exist, we do nothing */ if(parent_vertex_ptr!=NULL) { /* Get pointer to edge between the two vertices */ for(edge_ptr=EMANATING(vertex_ptr); edge_ptr!=NULL && VERTEX_PTR(edge_ptr)!=parent_vertex_ptr; edge_ptr=NEXT(edge_ptr)); if(edge_ptr==NULL) ens_error("Error:: couldn't find pointer to edge in reset_buses_branches\n"); /* Reset intensity */ EDGE_INTENSITY(edge_ptr) = Complex(0.,0.); } return(OK); } /*----------------------------------------------------------------------------------------------------*/ status init_bus_branch(graph_vertex vertex_ptr, graph_vertex parent_vertex_ptr, generic_ptr* input, generic_ptr* output) /* * It initializes the bus and preceding branch. * - The bus state is set to CONNECTED * - The bus voltage is set to 1+j0 pu * - The branch intensity is set to 0+j0 pu */ { graph_edge edge_ptr; VERTEX_STATE(vertex_ptr) = CONNECTED; VERTEX_VOLTAGE(vertex_ptr) = Complex(1.,0.); /* If the parent vertex does not exist, we do nothing */ if(parent_vertex_ptr!=NULL) { /* Get pointer to edge between the two vertices */ for(edge_ptr=EMANATING(vertex_ptr); edge_ptr!=NULL && VERTEX_PTR(edge_ptr)!=parent_vertex_ptr; edge_ptr=NEXT(edge_ptr)); if(edge_ptr==NULL) ens_error("Error:: couldn't find pointer to edge in reset_buses_branches\n"); /* Reset intensity */ EDGE_INTENSITY(edge_ptr) = Complex(0.,0.); } return(OK); } /****************************************************************************************************** * * Section 4 - Network manoeuvers * ******************************************************************************************************/ void toggle_branch_state_with_ptr(graph_edge edge_ptr) /* * Toggles the state of a branch (from its edge pointer) */ { switch(EDGE_STATE(edge_ptr)) { case OPEN: EDGE_STATE(edge_ptr) = CLOSED; EDGE_INTENSITY(edge_ptr) = Complex(0.,0.); break; case CLOSED: EDGE_STATE(edge_ptr)=OPEN; break; default: ens_error("Error in toggle_branch_state:: edge state can only be OPEN or CLOSED"); } } /*----------------------------------------------------------------------------------------------------*/ void toggle_branch_state_with_id(graph N, edge_id_type branch_id) /* * Toggles the state of a branch (from its edge id) */ { ens_error("Not implemented\n"); } /*----------------------------------------------------------------------------------------------------*/ void close_branch_with_id(graph N, edge_id_type branch_id) /* * Closes a branch (from its edge id) */ { set_state_branch_with_id(N,branch_id,CLOSED); } /*----------------------------------------------------------------------------------------------------*/ void open_branch_with_id(graph N, edge_id_type branch_id) /* * Opens a branch (from its edge id) */ { set_state_branch_with_id(N,branch_id,OPEN); } /*----------------------------------------------------------------------------------------------------*/ void close_branch_with_ptr(graph_edge edge_ptr) /* * Closes a branch (from its edge pointer) */ { ens_error("Not implemented\n"); } /*----------------------------------------------------------------------------------------------------*/ void open_branch_with_ptr(graph_edge edge_ptr) /* * Opens a branch from its edge pointer */ { ens_error("Not implemented\n"); } /*----------------------------------------------------------------------------------------------------*/ void set_state_branch_with_id(graph N, edge_id_type branch_id, int state) /* * Sets a branch state to "state" (from its edge id) */ { graph_edge edge_ptr; edge_ptr = find_branch_with_id(N, branch_id); if(edge_ptr!=NULL) { EDGE_STATE(edge_ptr)=state; if(state==OPEN) EDGE_INTENSITY(edge_ptr) = Complex(0.,0.); } else { ens_error("Error in set_state_branch_with_id:: branch not found"); } } /*----------------------------------------------------------------------------------------------------*/ void close_all_branches(graph N) /* * Closes all branches of the network N */ { graph_vertex v; graph_edge e; for(v=N;v!=NULL;v=NEXT(v)) for(e=EMANATING(v);e!=NULL;e=NEXT(e)) EDGE_STATE(e)=CLOSED; } /*----------------------------------------------------------------------------------------------------*/ void set_nominal_state(graph N) /* * Sets all switches of network N to their nominal state */ { graph_vertex v; graph_edge e; for(v=N;v!=NULL;v=NEXT(v)) for(e=EMANATING(v);e!=NULL;e=NEXT(e)) EDGE_STATE(e)=EDGE_STATE_NOMINAL(e); } /****************************************************************************************************** * * Section 5 - Miscellaneous * ******************************************************************************************************/ boolean branch_state_is_closed(graph_edge edge_ptr, graph_vertex vertex_ptr) /* * Tells whether the branch pointed to by edge_ptr is closed. */ { return(EDGE_STATE(edge_ptr)==CLOSED); } /*----------------------------------------------------------------------------------------------------*/ boolean branch_nominal_state_is_closed(graph_edge edge_ptr, graph_vertex vertex_ptr) /* * Tells whether the branch pointed to by edge_ptr is closed. */ { return(EDGE_STATE_NOMINAL(edge_ptr)==CLOSED); } /*----------------------------------------------------------------------------------------------------*/ void ens_error(char* error) /* * Reports an error */ { printf("%s\n", error); printf("Exiting...\n"); exit(1); } /*----------------------------------------------------------------------------------------------------*/ static boolean is_one_of(unsigned int num, unsigned int* numbers, int nnumbers) /* * Returns TRUE if the number 'num' is one the elements of the vector numbers[0..nnumbers]. */ { int i; for(i=0; iTINY_THRESH && fabs(z.i)>TINY_THRESH); }