network/node.cpp

// #define HYBRID // enables hybrid solver that can be used to perform state estimation (DPC 10/07: this is not fully functional)

#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <math.h>
#include "network.h"



// node CLASS FUNCTIONS

CLASS* node::oclass = NULL;
CLASS* node::pclass = NULL;
node *node::defaults = NULL;
//CLASS *node_class = NULL;

node::node(MODULE *mod) 
{

    // first time init
    if (oclass==NULL)
    {
        // register the class definition
        node_class = oclass = gl_register_class(mod,"node",PC_BOTTOMUP);
        if (oclass==NULL)
            GL_THROW("unable to register object class implemented by %s",__FILE__);

        // publish the class properties
        if (gl_publish_variable(oclass,
            PT_complex, "V", PADDR(V),
            PT_complex, "S", PADDR(S),
            PT_double, "G", PADDR(G),
            PT_double, "B", PADDR(B),
            PT_double, "Qmax_MVAR", PADDR(Qmax_MVAR),
            PT_double, "Qmin_MVAR", PADDR(Qmin_MVAR),
            PT_enumeration,"type",PADDR(type),
                PT_KEYWORD,"PQ",PQ,
                PT_KEYWORD,"PQV",PQV,
                PT_KEYWORD,"PV",PV,
                PT_KEYWORD,"SWING",SWING,
            PT_int16, "bus_id", PADDR(bus_id),
            PT_char32, "busname", PADDR(name),
            PT_int16, "flow_area_num", PADDR(flow_area_num),
            PT_double, "base_kV", PADDR(base_kV),
            PT_complex, "Vobs", PADDR(Vobs),
            PT_double, "Vstdev", PADDR(Vstdev),
            NULL)<1) GL_THROW("unable to publish properties in %s",__FILE__);

        // setup the default values
        defaults = this;
        type = PQ;
        V = complex(1,0,A);
        S = complex(0,0,J);
        G = 0;
        B = 0;
        Qmax_MVAR = 0.0;
        Qmin_MVAR = 0.0;
        Ys = 0.0;
        YVs = 0.0;
        type = PQ;
        name[0] = '\0';
        bus_id = 0;
        flow_area_num=1;
        loss_zone_num=1;
        base_kV=0;
        desired_kV=0;
        remote_bus_id=NULL;

        Vobs = complex(0,0,A); /* default observation is 0+0j */
        Vstdev = 0.0; /* default is no observation */
        r2 = 0.0;
        n_obs = 0;
        Sr2 = 0.0;
        n_inj = 0;
    }
}

node::~node()
{
}

int node::create() 
{
/*
    V = complex(1,0,A);
    S = complex(0,0,J);
    G = 0;
    B = 0;
    Qmax_MVAR = 0.0;
    Qmin_MVAR = 0.0;
    Ys = 0.0;
    YVs = 0.0;
    type = PQ;
    name[0] = '\0';
    bus_id = 0;
    flow_area_num=1;
    loss_zone_num=1;
    base_kV=0;
    desired_kV=0;
    remote_bus_id=NULL;

    Vobs = complex(0,0,A); // default observation is 0+0j 
    Vstdev = 0.0; // default is no observation 
    r2 = 0.0;
    n_obs = 0;
    Sr2 = 0.0;
    n_inj = 0;
*/
    memcpy(this,defaults,sizeof(*this));
    return 1;
}

int node::init(OBJECT *parent)
{
    // check that parent is swing bus
    OBJECT *hdr = OBJECTHDR(this);
    node *swing = parent?OBJECTDATA(parent,node):this;
    OBJECT *swing_hdr = OBJECTHDR(swing);
    if (swing_hdr->oclass!=hdr->oclass || swing->type!=SWING)
    {
        gl_error("node %s:%d parent is not a swing bus",hdr->oclass->name,hdr->id);
        return 0;
    }
#ifdef HYBRID
    // add observation residual
    if (Vstdev>0)
        swing->add_obs_residual(this);

    // add injection error
    swing->add_inj_residual(this);
#endif

    YVs=complex(0,0);
    Ys=complex(0,0);
    return 1;
}

#ifdef HYBRID
void node::add_inj_residual(node *pNode)
{
    if (pNode!=this)
    {
        pNode->Sr2 = (~(~pNode->V*((pNode->Ys + complex(pNode->G,pNode->B))*pNode->V + pNode->YVs)) - pNode->S).Mag();
        Sr2 += pNode->Sr2*pNode->Sr2;
        n_inj++;
    }
}

void node::del_inj_residual(node *pNode)
{
    if (pNode!=this)
    {
        Sr2 -= pNode->Sr2*pNode->Sr2;
        n_inj--;
    }
}

double node::get_inj_residual(void) const
{
    if (Sr2>0 && n_inj>0)
        return sqrt(Sr2/n_inj);
    else
        return Sr2;
}

void node::add_obs_residual(node *pNode)
{
    double r = (pNode->V - pNode->Vobs).Mag() / pNode->Vstdev;
    r*=r;
    if (pNode!=this)
    {
        pNode->r2 = r;
        pNode->n_obs = 1;
    }
    r2 += r;
    n_obs++;
}

void node::del_obs_residual(node *pNode)
{
    if (pNode!=this)
        r2 -= pNode->r2;
    else
    {
        double r = (pNode->V - pNode->Vobs).Mag() / pNode->Vstdev;
        r2 -= r*r;
    }
    n_obs--;
}

double node::get_obs_residual(void) const
{
    if (r2>0 && n_obs>0)
        return sqrt(r2/n_obs);
    else
        return r2;
}

double node::get_obs_probability(void) const
{
    if (r2<0)
        return 1;
    double pr = exp(-0.5*r2); 
    if (pr>1)
        gl_warning("node:%d observation probability exceeds 1!", OBJECTHDR(this)->id);
    return pr;
}
#endif

TIMESTAMP node::sync(TIMESTAMP t0) 
{
    OBJECT *hdr = OBJECTHDR(this);
    node *swing = hdr->parent?OBJECTDATA(hdr->parent,node):this;
    complex dV(0.0);
    complex YY = Ys + complex(G,B);
#ifdef HYBRID
    swing->del_inj_residual(this);
#endif
    if (!YY.IsZero() || type==SWING)
    {
        switch (type) {
        case PV:
            S.Im() = -((~V*(YY*V+YVs)).Im());
            if (Qmin_MVAR<Qmax_MVAR && S.Im()<Qmin_MVAR) 
                S.Im() = Qmin_MVAR;
            else if (Qmax_MVAR>Qmin_MVAR && S.Im()>Qmax_MVAR) 
                S.Im() = Qmax_MVAR;
            else
            {
                complex Vnew = (~S/~V - YVs) / YY;
                Vnew.SetPolar(V.Mag(),Vnew.Arg());
#ifdef HYBRID
                if (Vstdev>0)
                {
                    double pr = swing->get_obs_probability();
                    swing->del_obs_residual(this);
                    dV = Vobs*(1-pr) + (Vnew)*pr - V;
                    V += dV;
                    swing->add_obs_residual(this);
                }
                else
#endif
                {
                    dV = Vnew - V;
                    V = Vnew;
                }
                break;
            }
            /* continue with PQ solution */
        case PQ:
            if (!V.IsZero())
            {
                complex Vnew = (~S/~V - YVs) / YY;
#ifdef HYBRID
                if (Vstdev>0) // need to consider observation
                {
                    double pr = swing->get_obs_probability();
                    swing->del_obs_residual(this);
                    dV = Vobs*(1-pr) + (Vnew)*pr - V;
                    V += dV;
                    swing->add_obs_residual(this);
                }
                else // no observation 
#endif
                {
                    dV = (Vnew - V)*acceleration_factor;
                    V += dV;
                }
                V.Notation() = A;
            }
            break;
        case SWING:
            S = ~(~V*(YY*V + YVs));
            S.Notation() = J;
            break;
        default:
            /* unknown type fails */
            gl_error("invalid bus type");
            return TS_ZERO;
        }
    }
#ifdef HYBRID
    swing->add_inj_residual(this);
#endif

#ifdef _DEBUG
    // node debugging
    if (debug_node>0)
    {
        OBJECT* obj = OBJECTHDR(this);
        static int first=-1;
        if (first==-1) first = obj->id;
        if (obj->id==first)
        {
            printf("\n");
            printf("Node           Type  V                 Vobs              Stdev    Power             G        B        dV       Pr{Vobs} r2     Sr2\n");
            printf("============== ===== ================= ================= ======== ================= ======== ======== ======== ======== ====== ======\n");
        }
        if (((debug_node&1)==1 && dV.Mag()>convergence_limit )  // only on dV
#ifdef HYBRID
            || ((debug_node&2)==2 && Vstdev>0 ) // only on observation
            || ((debug_node&4)==4 && get_inj_residual()>0.001)// non-zero power residual
#endif
            )
        {
            printf("%2d (%-9.9s) %5s %+8.4f%+8.3fd %+8.4f%+8.3fd %8.5f %+8.4f%+8.4fj %+8.5f ", 
                obj->id, name, type==SWING?"SWING":(type==PQ?"PQ   ":"PV"),
                V.Mag(),V.Arg()*180/3.1416, 
                Vobs.Mag(), Vobs.Arg()*180/3.1416, Vstdev,
                S.Re(), S.Im(), 
                G, B,
                dV.Mag());
#ifdef HYBRID
            printf("%+8.5f %8.5f ", get_obs_probability(), r2);
            if (Vstdev>0)
                printf("%8.5f %6.3f %6.3f\n", get_obs_probability(), r2, get_inj_residual());
            else
                printf("   --      --   %6.3f\n",get_inj_residual());
#else
            printf("\n");
#endif
        }
    }
#endif // _DEBUG

    Ys = 0;
    YVs = 0;
    if (dV.Mag()>convergence_limit)
        return t0; /* did not converge, stop clock */
    else
        return TS_NEVER; /* converged, no further updates needed */
}

// IMPLEMENTATION OF CORE LINKAGE: node

EXPORT int create_node(OBJECT **obj, OBJECT *parent)
{
    *obj = gl_create_object(node_class,sizeof(node));
    if (*obj!=NULL)
    {
        last_node = *obj;
        node *my = OBJECTDATA(*obj,node);
        gl_set_parent(*obj,parent);
        my->create();
        return 1;
    }
    return 0;
}

EXPORT int init_node(OBJECT *obj)
{
    return OBJECTDATA(obj,node)->init(obj->parent);
}

EXPORT TIMESTAMP sync_node(OBJECT *obj, TIMESTAMP t0)
{
    TIMESTAMP t1 = OBJECTDATA(obj,node)->sync(t0);
    obj->clock = t0;
    return t1;
}

---

GridLAB-D^TM Version 1.0

An open-source project initiated by the US Department of Energy