GridLAB-D: powerflow/powerflow.cpp Source File

00001 
00025 #include <stdlib.h>
00026 #include <stdio.h>
00027 #include <errno.h>
00028 #include <math.h>
00029 
00031 double primary_voltage_ratio = 60.0;
00032 
00033 #include "powerflow.h"
00034 #include "node.h"
00035 #include "link.h"
00036 
00037 extern "C" {
00038 #include "output.h"
00039 }
00040 void print_matrix(complex mat[3][3])
00041 {
00042     for (int i = 0; i < 3; i++) {
00043         gl_testmsg("  %10.6f%+0.6f%c  %10.6f%+0.6f%c  %10.6f%+0.6f%c",
00044             mat[i][0].Re(), mat[i][0].Im(), mat[i][0].Notation(),
00045             mat[i][1].Re(), mat[i][1].Im(), mat[i][1].Notation(),
00046             mat[i][2].Re(), mat[i][2].Im(), mat[i][2].Notation());
00047     }
00048     gl_testmsg("\n");
00049 }
00050 
00051 
00052 EXPORT int kmldump(FILE *fp, OBJECT *obj)
00053 {
00054     if (obj==NULL) /* dump document styles */
00055     {
00056         /* line styles */
00057         fprintf(fp, "<Style id=\"overhead_line\">\n"
00058                     " <LineStyle>\n"
00059                     "  <color>7f00ffff</color>\n"
00060                     "  <width>4</width>\n"
00061                     " </LineStyle>\n"
00062                     " <PolyStyle>\n"
00063                     "  <color>7f00ff00</color>\n"
00064                     " </PolyStyle>\n"
00065                     "</Style>\n");
00066         fprintf(fp, "<Style id=\"underground_line\">\n"
00067                     " <LineStyle>\n"
00068                     "  <color>3f00ffff</color>\n"
00069                     "  <width>4</width>\n"
00070                     " </LineStyle>\n"
00071                     " <PolyStyle>\n"
00072                     "  <color>3f00ff00</color>\n"
00073                     " </PolyStyle>\n"
00074                     "</Style>\n");
00075     }
00077     else if (gl_object_isa(obj,"node"))
00078     {
00079         node *pNode = OBJECTDATA(obj,node);
00080         if (isnan(obj->latitude) || isnan(obj->longitude))
00081             return 0;
00082         fprintf(fp,"<Placemark>\n");
00083         if (obj->name)
00084             fprintf(fp,"<name>%s</name>\n", obj->name, obj->oclass->name, obj->id);
00085         else
00086             fprintf(fp,"<name>%s %d</name>\n", obj->oclass->name, obj->id);
00087         fprintf(fp,"<description>\n");
00088         fprintf(fp,"<![CDATA[\n");
00089         fprintf(fp,"<TABLE>\n");
00090         fprintf(fp,"<TR><TD WIDTH=\"25%\">%s&nbsp;%d<HR></TD><TH WIDTH=\"25%\" ALIGN=CENTER>Phase A<HR></TH><TH WIDTH=\"25%\" ALIGN=CENTER>Phase B<HR></TH><TH WIDTH=\"25%\" ALIGN=CENTER>Phase C<HR></TH></TR>\n", obj->oclass->name, obj->id);
00091 
00092         // voltages
00093         fprintf(fp,"<TR><TH ALIGN=LEFT>Voltage</TH>");
00094         double vscale = primary_voltage_ratio*sqrt((double) 3.0)/(double) 1000.0;
00095         if (pNode->has_phase(PHASE_A))
00096             fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;kV&nbsp;&nbsp;<BR>%.3f&nbsp;deg&nbsp;</TD>",
00097                 pNode->phaseA_V.Mag()*vscale,pNode->phaseA_V.Arg()*180/3.1416);
00098         else
00099             fprintf(fp,"<TD></TD>");
00100         if (pNode->has_phase(PHASE_B))
00101             fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;kV&nbsp;&nbsp;<BR>%.3f&nbsp;deg&nbsp;</TD>",
00102                 pNode->phaseB_V.Mag()*vscale,pNode->phaseB_V.Arg()*180/3.1416);
00103         else
00104             fprintf(fp,"<TD></TD>");
00105         if (pNode->has_phase(PHASE_C))
00106             fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;kV&nbsp;&nbsp;<BR>%.3f&nbsp;deg&nbsp;</TD>",
00107                 pNode->phaseC_V.Mag()*vscale,pNode->phaseC_V.Arg()*180/3.1416);
00108         else
00109             fprintf(fp,"<TD></TD>");
00110 
00111         // supply
00113 
00114         // demand
00115         if (gl_object_isa(obj,"load"))
00116         {
00117             load *pLoad = OBJECTDATA(obj,load);
00118             fprintf(fp,"<TR><TH ALIGN=LEFT>Load</TH>");
00119             if (pNode->has_phase(PHASE_A))
00120             {
00121                 complex load_A = ~pNode->phaseA_V*pLoad->phaseA_I + pLoad->phaseA_VA;
00122                 fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;kW&nbsp;&nbsp;<BR>%.3f&nbsp;kVAR</TD>",
00123                     load_A.Re(),load_A.Im());
00124             }
00125             else
00126                 fprintf(fp,"<TD></TD>");
00127             if (pNode->has_phase(PHASE_B))
00128             {
00129                 complex load_B = ~pNode->phaseB_V*pLoad->phaseB_I + pLoad->phaseB_VA;
00130                 fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;kW&nbsp;&nbsp;<BR>%.3f&nbsp;kVAR</TD>",
00131                     load_B.Re(),load_B.Im());
00132             }
00133             else
00134                 fprintf(fp,"<TD></TD>");
00135             if (pNode->has_phase(PHASE_C))
00136             {
00137                 complex load_C = ~pNode->phaseC_V*pLoad->phaseC_I + pLoad->phaseC_VA;
00138                 fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;kW&nbsp;&nbsp;<BR>%.3f&nbsp;kVAR</TD>",
00139                     load_C.Re(),load_C.Im());
00140             }
00141             else
00142                 fprintf(fp,"<TD></TD>");
00143         }
00144         fprintf(fp,"</TR>\n");
00145         fprintf(fp,"</TABLE>\n");
00146         fprintf(fp,"]]>\n");
00147         fprintf(fp,"</description>\n");
00148         fprintf(fp,"<Point>\n");
00149         fprintf(fp,"<coordinates>%f,%f</coordinates>\n",obj->longitude,obj->latitude);
00150         fprintf(fp,"</Point>\n");
00151         fprintf(fp,"</Placemark>\n");
00152         return 0;
00153     }
00154     else if (gl_object_isa(obj,"link"))
00155     {
00156         link *pLink = OBJECTDATA(obj,link);
00157         OBJECT *from = pLink->from;
00158         OBJECT *to = pLink->to;
00159         if (isnan(from->latitude) || isnan(to->latitude) || isnan(from->longitude) || isnan(to->longitude))
00160             return 0;
00161         fprintf(fp,"    <Placemark>\n");
00162         if (obj->name)
00163             fprintf(fp,"      <name>%s</name>\n", obj->name);
00164         else
00165             fprintf(fp,"      <name>%s ==> %s</name>\n", from->name?from->name:"unnamed", to->name?to->name:"unnamed");
00166         fprintf(fp,"      <description>\n");
00167         fprintf(fp,"        <![CDATA[\n");
00168         fprintf(fp,"          <TABLE><TR>\n");
00169         fprintf(fp,"<TR><TD WIDTH=\"25%\">%s&nbsp;%d<HR></TD><TH WIDTH=\"25%\" ALIGN=CENTER>Phase A<HR></TH><TH WIDTH=\"25%\" ALIGN=CENTER>Phase B<HR></TH><TH WIDTH=\"25%\" ALIGN=CENTER>Phase C<HR></TH></TR>\n", obj->oclass->name, obj->id);
00170 
00171         // values
00172         node *pFrom = OBJECTDATA(from,node);
00173         node *pTo = OBJECTDATA(to,node);
00174         double vscale = primary_voltage_ratio*sqrt((double) 3.0)/(double) 1000.0;
00175         complex loss[3];
00176         complex flow[3];
00177         complex current[3];
00178         complex in[3] = {pLink->i_abc_in[0], pLink->i_abc_in[1], pLink->i_abc_in[2]};
00179         complex out[3] = {pLink->i_abc_out[0], pLink->i_abc_out[1], pLink->i_abc_out[2]};
00180         complex Vfrom[3] = {pFrom->phaseA_V, pFrom->phaseB_V, pFrom->phaseC_V};
00181         complex Vto[3] = {pTo->phaseA_V,pTo->phaseB_V,pTo->phaseC_V};
00182         int phase[3] = {pLink->has_phase(PHASE_A),pLink->has_phase(PHASE_B),pLink->has_phase(PHASE_C)};
00183         int i;
00184         for (i=0; i<3; i++)
00185         {
00186             if (phase[i])
00187             {
00188                 if (Vfrom[i].Re() > Vto[i].Re())
00189                 {
00190                     flow[i] = out[i]*Vto[i]*vscale;
00191                     loss[i] = in[i]*Vfrom[i]*vscale - flow[i];
00192                     current[i] = out[i];
00193                 }
00194                 else
00195                 {
00196                     flow[i] = in[i]*Vfrom[i]*vscale;
00197                     loss[i] = out[i]*Vto[i]*vscale - flow[i];
00198                     current[i] = in[i];
00199                 }
00200             }
00201         }
00202 
00203         // flow
00204         fprintf(fp,"<TR><TH ALIGN=LEFT>Flow</TH>");
00205         for (i=0; i<3; i++)
00206         {
00207             if (phase[i])
00208                 fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;&nbsp;kW&nbsp;&nbsp;<BR>%.3f&nbsp;&nbsp;kVAR</TD>\n",
00209                     flow[i].Re(), flow[i].Im());
00210             else
00211                 fprintf(fp,"<TD></TD>\n");
00212         }
00213         fprintf(fp,"</TR>");
00214 
00215         // current
00216         fprintf(fp,"<TR><TH ALIGN=LEFT>Current</TH>");
00217         for (i=0; i<3; i++)
00218         {
00219             if (phase[i])
00220                 fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.3f&nbsp;&nbsp;Amps</TD>\n",
00221                     current[i].Mag());
00222             else
00223                 fprintf(fp,"<TD></TD>\n");
00224         }
00225         fprintf(fp,"</TR>");
00226 
00227         // loss
00228         fprintf(fp,"<TR><TH ALIGN=LEFT>Loss</TH>");
00229         for (i=0; i<3; i++)
00230         {
00231             if (phase[i])
00232                 fprintf(fp,"<TD ALIGN=RIGHT STYLE=\"font-family:courier;\">%.2f&nbsp;&nbsp;&nbsp;%%P&nbsp;&nbsp;<BR>%.2f&nbsp;&nbsp;&nbsp;%%Q&nbsp;&nbsp;</TD>\n",
00233                     loss[i].Re()/flow[i].Re()*100,loss[i].Im()/flow[i].Im()*100);
00234             else
00235                 fprintf(fp,"<TD></TD>\n");
00236         }
00237         fprintf(fp,"</TR>");
00238         fprintf(fp,"</TABLE>\n");
00239         fprintf(fp,"        ]]>\n");
00240         fprintf(fp,"      </description>\n");
00241         fprintf(fp,"      <styleUrl>#%s</styleUrl>>\n",obj->oclass->name);
00242         fprintf(fp,"      <coordinates>%f,%f</coordinates>\n",
00243             (from->longitude+to->longitude)/2,(from->latitude+to->latitude)/2);
00244         fprintf(fp,"      <LineString>\n");         
00245         fprintf(fp,"        <extrude>0</extrude>\n");
00246         fprintf(fp,"        <tessellate>0</tessellate>\n");
00247         fprintf(fp,"        <altitudeMode>relative</altitudeMode>\n");
00248         fprintf(fp,"        <coordinates>%f,%f,50 %f,%f,50</coordinates>\n",
00249             from->longitude,from->latitude,to->longitude,to->latitude);
00250         fprintf(fp,"      </LineString>\n");            
00251         fprintf(fp,"    </Placemark>\n");
00252         return 0;
00253     }
00254     return 0; /* 1 means output default if it wasn't handled */
00255 }
00256