powerflow/regulator.cpp

// $Id: regulator.cpp,v 1.17 2008/02/13 22:48:04 natet Exp $
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <math.h>
#include <iostream>
using namespace std;

#include "regulator.h"
#include "node.h"

CLASS* regulator::oclass = NULL;
CLASS* regulator::pclass = NULL;
regulator *regulator::defaults = NULL;

regulator::regulator(MODULE *mod) : link(mod)
{
    if (oclass==NULL)
    {
        // link to parent class (used by isa)
        pclass = link::oclass;

        // register the class definition
        oclass = gl_register_class(mod,"regulator",PC_BOTTOMUP|PC_POSTTOPDOWN);
        if (oclass==NULL)
            GL_THROW("unable to register object class implemented by %s",__FILE__);

        // publish the class properties
        if (gl_publish_variable(oclass,
            PT_object,"configuration",PADDR(configuration),
            PT_set, "phases", PADDR(phases),
               PT_KEYWORD, "A",PHASE_A,
               PT_KEYWORD, "B",PHASE_B,
               PT_KEYWORD, "C",PHASE_C,
               PT_KEYWORD, "N",PHASE_N,
               PT_KEYWORD, "S",PHASE_S,
            PT_object, "from",PADDR(from),
            PT_object, "to", PADDR(to),
            NULL)<1) GL_THROW("unable to publish properties in %s",__FILE__);

        // setup the default values
        defaults = this;
        memset(this,0,sizeof(regulator));
    }
}

int regulator::isa(char *classname)
{
    return strcmp(classname,"regulator")==0 || link::isa(classname);
}

int regulator::create()
{
    int result = link::create();
    configuration = NULL;
    return result;
}

int regulator::init()
{
    int result = link::init();

    if (!configuration)
        throw "no regulator configuration specified.";

    if (!gl_object_isa(configuration, "regulator_configuration"))
        throw "invalid regulator configuration";

    regulator_configuration *config = OBJECTDATA(configuration, regulator_configuration);
    node *volt_node = OBJECTDATA(to, node);
    //tap_pos = {0, 0, 0};     // tap position
    complex D_mat[3][3] = {{complex(1,0),complex(-1,0),complex(0,0)},
                           {complex(0,0), complex(1,0),complex(-1,0)},
                           {complex(-1,0), complex(0,0), complex(1,0)}};   // D_mat - 3x3 matrix, 'D' matrix
    complex W_mat[3][3] = {{complex(2,0),complex(1,0),complex(0,0)},
                           {complex(0,0),complex(2,0),complex(1,0)},
                           {complex(1,0),complex(0,0),complex(2,0)}};
    
    multiply(1.0/3.0,W_mat,W_mat);
    complex volt[3] = {0.0, 0.0, 0.0};
    double tapChangePer = config->regulation / (double) config->raise_taps;
    double VtapChange = (config->V_high / ((double) config->PT_ratio * sqrt(3.0))) * tapChangePer;
    double Vlow = config->band_center - config->band_width / 2.0;
    double Vhigh = config->band_center + config->band_width / 2.0;
    complex V2[3], Vcomp[3];

    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            a_mat[i][j] = b_mat[i][j] = c_mat[i][j] = d_mat[i][j] =
                    A_mat[i][j] = B_mat[i][j] = 0.0;
        }
        //config->tap_pos[i] = 0;
        i_abc_in[i] = i_n_in[i] = i_abc_out[i] = i_n_out[i] = 1.0;
    }

    if (volt_node) {
        volt[0] = volt_node->phaseA_V;
        volt[1] = volt_node->phaseB_V;
        volt[2] = volt_node->phaseC_V;
    }
    else
        volt[0] = volt[1] = volt[2] = 0.0;

    for (int i = 0; i < 3; i++) {
        V2[i] = volt[i] / ((double) config->PT_ratio * sqrt(3.0));
        Vcomp[i] = V2[i] - (i_abc_out[i] / (double) config->CT_ratio) * complex(config->ldc_R_V, config->ldc_X_V);

        /*if (Vcomp[i].Re() < Vlow)
            config->tap_pos[i] += (int16) floor(fabs(Vlow - V2[i].Re()) / VtapChange);
        else if (Vcomp[i].Re() > Vhigh)
            config->tap_pos[i] += (int16) floor(fabs(Vhigh - V2[i].Re()) / VtapChange);
        */
        a_mat[i][i] = 1.0 - config->tap_pos[i] * tapChangePer;
    }

    complex tmp_mat[3][3] = {{complex(1,0)/a_mat[0][0],complex(0,0),complex(0,0)},
                             {complex(0,0), complex(1,0)/a_mat[1][1],complex(0,0)},
                             {complex(-1,0)/a_mat[0][0],complex(-1,0)/a_mat[1][1],complex(0,0)}};
    complex tmp_mat1[3][3];

    switch ((enum connect_type) config->connect_type) {
        case WYE_WYE:
            for (int i = 0; i < 3; i++)
                d_mat[i][i] = complex(1.0,0) / a_mat[i][i]; // forced to use complex math.  I hope this works!
                inverse(a_mat,A_mat);
            break;
        case OPEN_DELTA_ABBC:
            //for (int i = 0; i < 3; i++)
                //d_mat[i][i] = a_mat[i][i];
            d_mat[0][0] = complex(1,0) / a_mat[0][0];
            d_mat[1][0] = complex(-1,0) / a_mat[0][0];
            d_mat[1][2] = complex(-1,0) / a_mat[1][1];
            d_mat[2][2] = complex(1,0) / a_mat[1][1];

            a_mat[2][0] = -a_mat[0][0];
            a_mat[2][1] = -a_mat[1][1];
            a_mat[2][2] = 0;

            // Fix tmp matrix
            //gl_testmsg("regulator: W matrix\n");
            //print_matrix(W_mat);
            //gl_testmsg("regulator: D matrix\n");
            //print_matrix(D_mat);
            //print_matrix(tmp_mat);
            multiply(W_mat,tmp_mat,tmp_mat1);
            //print_matrix(tmp_mat1);
            multiply(tmp_mat1,D_mat,A_mat);
            break;
        case OPEN_DELTA_BCAC:
            break;
        case OPEN_DELTA_CABA:
            break;
        case CLOSED_DELTA:

            break;
        default:
            throw "unknown regulator connect type";
            break;
    }
#ifdef _TESTING
    gl_testmsg("regulator: a matrix\n");
    print_matrix(a_mat);
    gl_testmsg("regulator: A matrix\n");
    print_matrix(A_mat);
    gl_testmsg("regulator: b matrix\n");
    print_matrix(b_mat);
    gl_testmsg("regulator: B matrix\n");
    print_matrix(B_mat);
    gl_testmsg("regulator: d matrix\n");
    print_matrix(d_mat);
#endif
    return result;
}
// IMPLEMENTATION OF CORE LINKAGE: regulator

EXPORT int create_regulator(OBJECT **obj, OBJECT *parent)
{
    try
    {
        *obj = gl_create_object(regulator::oclass,sizeof(regulator));
        if (*obj!=NULL)
        {
            regulator *my = OBJECTDATA(*obj,regulator);
            gl_set_parent(*obj,parent);
            return my->create();
        }
    }
    catch (char *msg)
    {
        gl_error("create_regulator: %s", msg);
    }
    return 1;
}

EXPORT int init_regulator(OBJECT *obj)
{
    regulator *my = OBJECTDATA(obj,regulator);
    try {
        return my->init();
    }
    catch (char *msg)
    {
        GL_THROW("%s (regulator:%d): %s", my->get_name(), my->get_id(), msg);
        return 0; 
    }
}

EXPORT TIMESTAMP sync_regulator(OBJECT *obj, TIMESTAMP t0, PASSCONFIG pass)
{
    regulator *pObj = OBJECTDATA(obj,regulator);
    try {
        if (pass==PC_BOTTOMUP)
            return pObj->sync(t0);
        else
        {
            TIMESTAMP t1 = pObj->postsync(t0);
            obj->clock = t0;
            return t1;
        }
    } catch (const char *error) {
        GL_THROW("%s (%s:%d): %s", pObj->get_name(), pObj->get_id(), error);
        return 0; 
    } catch (...) {
        GL_THROW("%s (%s:%d): %s", pObj->get_name(), pObj->get_id(), "unknown exception");
        return 0;
    }
}

EXPORT int isa_regulator(OBJECT *obj, char *classname)
{
    return OBJECTDATA(obj,regulator)->isa(classname);
}

---

GridLAB-D^TM Version 1.0

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