commercial/office.cpp

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <math.h>
#include "gridlabd.h"
#include "solvers.h"
#include "office.h"
#include "../climate/climate.h"

/* the object class definition */
CLASS *office::oclass = NULL;

/* the static default object properties */
office *office::defaults = NULL;

/* specify passes that are needed */
static PASSCONFIG passconfig = PC_PRETOPDOWN|PC_BOTTOMUP;

/* specify pass that advances the clock  */
static PASSCONFIG clockpass = PC_BOTTOMUP;

/* Class registration is only called once to register the class with the core */
office::office(MODULE *module)
{
    if (oclass==NULL)
    {
        oclass = gl_register_class(module,"office",passconfig); 
        if (gl_publish_variable(oclass,
            PT_double, "occupancy", PADDR(zone.current.occupancy),
            PT_double, "floor_height[ft]", PADDR(zone.design.floor_height),
            PT_double, "exterior_ua[Btu/degF/h]", PADDR(zone.design.exterior_ua),
            PT_double, "interior_ua[Btu/degF/h]", PADDR(zone.design.interior_ua),
            PT_double, "interior_mass[Btu/degF]", PADDR(zone.design.interior_mass),
            PT_double, "floor_area[sf]", PADDR(zone.design.floor_area),
            PT_double, "glazing[sf]", PADDR(zone.design.window_area[0]), PT_SIZE, sizeof(zone.design.window_area)/sizeof(zone.design.window_area[0]),
            PT_double, "glazing.north[sf]", PADDR(zone.design.window_area[CP_N]),
            PT_double, "glazing.northeast[sf]", PADDR(zone.design.window_area[CP_NE]),
            PT_double, "glazing.east[sf]", PADDR(zone.design.window_area[CP_E]),
            PT_double, "glazing.southeast[sf]", PADDR(zone.design.window_area[CP_SE]),
            PT_double, "glazing.south[sf]", PADDR(zone.design.window_area[CP_S]),
            PT_double, "glazing.southwest[sf]", PADDR(zone.design.window_area[CP_SW]),
            PT_double, "glazing.west[sf]", PADDR(zone.design.window_area[CP_W]),
            PT_double, "glazing.northwest[sf]", PADDR(zone.design.window_area[CP_NW]),
            PT_double, "glazing.horizontal[sf]", PADDR(zone.design.window_area[CP_H]),
            PT_double, "glazing.coefficient[pu]", PADDR(zone.design.glazing_coeff),

            PT_double, "occupants", PADDR(zone.design.occupants),
            PT_char256, "schedule", PADDR(zone.design.schedule),

            PT_double, "air_temperature[degF]", PADDR(zone.current.air_temperature),
            PT_double, "mass_temperature[degF]", PADDR(zone.current.mass_temperature),
            PT_double, "temperature_change[degF/h]", PADDR(zone.current.temperature_change),
            PT_double, "Qh[Btu/h]", PADDR(Qh),
            PT_double, "Qs[Btu/h]", PADDR(Qs),
            PT_double, "Qi[Btu/h]", PADDR(Qi),
            PT_double, "Qz[Btu/h]", PADDR(Qz),

            /* HVAC loads */
            PT_enumeration, "hvac.mode", PADDR(zone.hvac.mode),
                PT_KEYWORD, "HEAT", HC_HEAT,
                PT_KEYWORD, "AUX", HC_AUX,
                PT_KEYWORD, "COOL", HC_COOL,
                PT_KEYWORD, "ECON", HC_ECON,
                PT_KEYWORD, "OFF", HC_OFF,
            PT_double, "hvac.cooling.balance_temperature[degF]", PADDR(zone.hvac.cooling.balance_temperature),
            PT_double, "hvac.cooling.capacity[Btu/h]", PADDR(zone.hvac.cooling.capacity),
            PT_double, "hvac.cooling.capacity_perF[Btu/degF/h]",PADDR(zone.hvac.cooling.capacity_perF),
            PT_double, "hvac.cooling.design_temperature[degF]", PADDR(zone.hvac.cooling.design_temperature),
            PT_double, "hvac.cooling.efficiency[pu]",PADDR(zone.hvac.cooling.efficiency),
            PT_double, "hvac.cooling.cop[pu]", PADDR(zone.hvac.cooling.cop),
            PT_double, "hvac.heating.balance_temperature[degF]",PADDR(zone.hvac.heating.balance_temperature),
            PT_double, "hvac.heating.capacity[Btu/h]",PADDR(zone.hvac.heating.capacity),
            PT_double, "hvac.heating.capacity_perF[Btu/degF/h]", PADDR(zone.hvac.heating.capacity_perF),
            PT_double, "hvac.heating.design_temperature[degF]", PADDR(zone.hvac.heating.design_temperature),
            PT_double, "hvac.heating.efficiency[pu]", PADDR(zone.hvac.heating.design_temperature),
            PT_double, "hvac.heating.cop[pu]", PADDR(zone.hvac.heating.cop),

            /* Lighting loads */
            PT_double, "lights.capacity[kW]", PADDR(zone.lights.capacity),
            PT_double, "lights.fraction[pu]", PADDR(zone.lights.fraction),

            /* Plug loads */
            PT_double, "plugs.capacity[kW]", PADDR(zone.plugs.capacity),
            PT_double, "plugs.fraction[pu]", PADDR(zone.plugs.fraction),

            /* End-use data */
            PT_complex, "demand[kW]", PADDR(zone.total.demand),
            PT_complex, "power[kW]", PADDR(zone.total.power),
            PT_complex, "energy", PADDR(zone.total.energy),
            PT_double, "power_factor", PADDR(zone.total.power_factor),

            PT_complex, "hvac.demand[kW]", PADDR(zone.hvac.enduse.demand),
            PT_complex, "hvac.power[kW]", PADDR(zone.hvac.enduse.power),
            PT_complex, "hvac.energy", PADDR(zone.hvac.enduse.energy),
            PT_double, "hvac.power_factor", PADDR(zone.hvac.enduse.power_factor),

            PT_complex, "lights.demand[kW]", PADDR(zone.lights.enduse.demand),
            PT_complex, "lights.power[kW]", PADDR(zone.lights.enduse.power),
            PT_complex, "lights.energy", PADDR(zone.lights.enduse.energy),
            PT_double, "lights.power_factor", PADDR(zone.lights.enduse.power_factor),

            PT_complex, "plugs.demand[kW]", PADDR(zone.plugs.enduse.demand),
            PT_complex, "plugs.power[kW]", PADDR(zone.plugs.enduse.power),
            PT_complex, "plugs.energy", PADDR(zone.plugs.enduse.energy),
            PT_double, "plugs.power_factor", PADDR(zone.plugs.enduse.power_factor),

            PT_double, "control.cooling_setpoint", PADDR(zone.control.cooling_setpoint),
            PT_double, "control.heating_setpoint", PADDR(zone.control.heating_setpoint),
            PT_double, "control.setpoint_deadband", PADDR(zone.control.setpoint_deadband),
            PT_double, "control.ventilation_fraction", PADDR(zone.control.ventilation_fraction),
            PT_double, "control.lighting_fraction", PADDR(zone.control.lighting_fraction),
            
            NULL)<1) throw("unable to publish properties in "__FILE__);
        defaults = this;
        memset(defaults,0,sizeof(office));

        /* set default power factors */
        zone.lights.enduse.power_factor = 1.0;
        zone.plugs.enduse.power_factor = 1.0;
        zone.hvac.enduse.power_factor = 1.0;

        /* set default climate to static values */
        static double Tout = 59, RHout=0.75, Solar[9]={1, 0.9,0.9, 0.5,0.5, 0,0, 0,1};
        zone.current.pTemperature = &Tout;
        zone.current.pHumidity = &RHout;
        zone.current.pSolar = Solar;

        /* set default thermal conditions */
        zone.current.air_temperature = Tout;
        zone.current.mass_temperature = Tout;

        /* set default control strategy */
        zone.control.heating_setpoint = 70;
        zone.control.cooling_setpoint = 75;
        zone.control.setpoint_deadband = 1;
        zone.control.ventilation_fraction = 0.2;
        zone.control.lighting_fraction = 0.5;
    }
}

/* Object creation is called once for each object that is created by the core */
int office::create(void) 
{
    memcpy(this,defaults,sizeof(*this));
    return 1; /* return 1 on success, 0 on failure */
}

/* Object initialization is called once after all object have been created */
int office::init(OBJECT *parent)
{
    update_control_setpoints();

    /* automatic sizing of equipment */
    if (zone.hvac.heating.capacity==0)
        zone.hvac.heating.capacity = zone.design.exterior_ua*(zone.control.heating_setpoint-zone.hvac.heating.design_temperature);
    if (zone.hvac.cooling.capacity==0)
        zone.hvac.cooling.capacity = -zone.design.exterior_ua*(zone.hvac.cooling.design_temperature-zone.control.cooling_setpoint) /* losses */
            - (zone.design.window_area[0]+zone.design.window_area[1]+zone.design.window_area[2]+zone.design.window_area[3]+zone.design.window_area[4]
                +zone.design.window_area[5]+zone.design.window_area[6]+zone.design.window_area[7]+zone.design.window_area[8])*100*3.412*zone.design.glazing_coeff /* solar */
            - (zone.lights.capacity + zone.plugs.capacity)*3.412;
    if (zone.hvac.cooling.cop==0)
        zone.hvac.cooling.cop=-3.5;
    if (zone.hvac.heating.cop==0)
        zone.hvac.heating.cop=1.5;

    OBJECT *hdr = OBJECTHDR(this);

    // link to climate data
    static FINDLIST *climates = gl_find_objects(FL_NEW,FT_CLASS,SAME,"climate",FT_END);
    if (climates==NULL)
        gl_warning("house: no climate data found, using static data");
    else if (climates->hit_count>1)
        gl_warning("house: %d climates found, using first one defined", climates->hit_count);
    if (climates->hit_count>0)
    {
        OBJECT *obj = gl_find_next(climates,NULL);
        if (obj->rank<=hdr->rank)
            gl_set_dependent(obj,hdr);
        zone.current.pTemperature = (double*)GETADDR(obj,gl_get_property(obj,"temperature"));
        zone.current.pHumidity = (double*)GETADDR(obj,gl_get_property(obj,"humidity"));
        zone.current.pSolar = (double*)GETADDR(obj,gl_get_property(obj,"solar_flux"));
    }

    struct {
        char *desc;
        bool test;
    } map[] = {
        {"floor height is not valid", zone.design.floor_height<=0},
        {"interior mass is not valid", zone.design.interior_mass<=0},
        {"interior UA is not valid", zone.design.interior_ua<=0},
        {"exterior UA is not valid", zone.design.exterior_ua<=0},
        {"floor area is not valid",zone.design.floor_area<=0},
        {"control setpoint deadpoint is invalid", zone.control.setpoint_deadband<=0},
        {"heating and cooling setpoints conflict",TheatOn>=TcoolOff},
        {"heating balance temperature is not greater than heating design temperature",zone.hvac.heating.design_temperature>=zone.hvac.heating.balance_temperature},
        {"cooling balance temperature is not less than cooling design temperature",zone.hvac.cooling.design_temperature<=zone.hvac.cooling.balance_temperature},
        {"cooling capacity is not negative", zone.hvac.cooling.capacity>=0},
        {"heating capacity is not positive", zone.hvac.heating.capacity<=0},
        {"cooling cop is not negative", zone.hvac.cooling.cop>=0},
        {"heating cop is not positive", zone.hvac.heating.cop<=0},
    };
    int i;
    for (i=0; i<sizeof(map)/sizeof(map[0]); i++)
    {
        if (map[i].test)
            throw map[i].desc;
    }
    return 1; /* return 1 on success, 0 on failure */
}

/* Sync is called when the clock needs to advance on the bottom-up pass */
TIMESTAMP office::presync(TIMESTAMP t0, TIMESTAMP t1) 
{
    /* reset the multizone heat transfer */
    Qz = 0;
    return TS_NEVER;
}

/* Sync is called when the clock needs to advance on the bottom-up pass */
TIMESTAMP office::sync(TIMESTAMP t0, TIMESTAMP t1) 
{
    /* local aliases */
    const double &Tout = (*(zone.current.pTemperature));
    const double &UA = (zone.design.exterior_ua);
    const double &Cm = (zone.design.interior_mass);
    const double &Hm = (zone.design.interior_ua);
    double &Ti = (zone.current.air_temperature);
    double &dTi = (zone.current.temperature_change);
    double &Tm = (zone.current.mass_temperature);
    HCMODE &mode = (zone.hvac.mode);

    /* advance the thermal state of the building */
    const double dt1 = t0>0 ? (double)(t1-t0)*TS_SECOND : 0;
    if (dt1>0)
    {
        const double Ca = 0.2402 * 0.0735 * zone.design.floor_height * zone.design.floor_area;
        const double dt = dt1/3600;

        /* update enduses and get internal heat gains */
        Qi = update_lighting(dt) + update_plugs(dt);

        /* compute heating/cooling effect */
        Qh = update_hvac(dt); 

        Qs = 0; 
        int i;
        for (i=0; i<9; i++)
            Qs += 3.412 * zone.design.window_area[i] * zone.current.pSolar[i]/10;
        Qs *= dt;
        if (Qs<0)
            throw "solar gain is negative?!?";

        /* compute thermal constants */
        const double c1 = -(UA + Hm)/Ca;
        const double c2 = Hm/Ca;
        const double c3 = (Qi + Qh + Qs + Qz + Tout*UA)/Ca;
        const double c4 = Hm/Cm;
        const double c5 = -c4;
        const double p1 = 1/c2;
        const double p2 = -(c5+c1)/c2;
        const double p3 = c1*c5/c2 - c4;
        const double p4 = -c3*c5/c2;

        Teq = p4/p3;
        const double ra = 2*p1;
        const double rb = -p2/ra;
        const double rc = sqrt(p2*p2-4*p1*p3)/ra;
        if (!isfinite(rc))
            throw "thermal solution does not have any roots";
        r1 = rb+rc;
        r2 = rb-rc;
        if (r1==r2)
            throw "thermal solution has only one root";

        /* compute initial condition constants */
        dTi = c2*Tm + c1*Ti + c3;
        k2 = (dTi-r1*Ti+r1*Teq)/(r2-r1);
        k1 = Ti-k2-Teq;

        /* compute the next temperature shift */
        const double ker1 = k1*exp(r1*dt);
        const double ker2 = k2*exp(r2*dt);
        const double Ts = ker1+ker2;
        if (!isfinite(Ts))
            throw("temperature balance equation failed");

        /* compute the new dependent values */
        Ti = Teq+Ts;
        dTi = r1*ker1+r2*ker2;
        Tm = (dTi-c1*Ti-c3)/c2;

        /* calculate the power consumption */
        zone.total.power = zone.lights.enduse.power + zone.plugs.enduse.power + zone.hvac.enduse.power;
        zone.total.energy = zone.total.power * dt;
    }

    /* determine the temperature of the next event */
#define TPREC 0.01
    unsigned int Nevents = 0;
    double dt2=(double)TS_NEVER;
    Tevent = Teq;
    if (Qh<0) /* cooling active */
        Tevent = TcoolOff;
    else if (Qh>0) /* heating active */
        Tevent = TheatOff;
    else if (Teq<=TheatOn+TPREC) /* heating deadband */
        Tevent = TheatOn;
    else if (Teq>=TcoolOn-TPREC) /* cooling deadband */
        Tevent = TcoolOn;
    else /* equilibrium */
        return TS_NEVER;

    /* solve for the time to the next event */
    Nevents = dual_decay_solve(&dt2,TPREC,0,dt2,k1,r1,k2,r2,Teq-Tevent);
    if (Nevents<0)
        throw("dual_decay_solve() failed");
    else if (Nevents==0)
        return TS_NEVER;
    else if (dt2==0) /* first solution is immediate, advance clock 1 tick to get next solution */
        dt2=1.0/3600/TS_SECOND; /* advance just 1 tick */
    if (dt2==TS_NEVER)
        return TS_NEVER;
    else
        return t1+(TIMESTAMP)(dt2*3600*TS_SECOND); /* return t2>t1 on success, t2=t1 for retry, t2<t1 on failure */
}

void office::update_control_setpoints()
{
    TcoolOn = zone.control.cooling_setpoint + zone.control.setpoint_deadband;
    TcoolOff = zone.control.cooling_setpoint - zone.control.setpoint_deadband;
    TheatOn = zone.control.heating_setpoint - zone.control.setpoint_deadband;
    TheatOff = zone.control.heating_setpoint + zone.control.setpoint_deadband;
}

double office::update_lighting(double dt)
{
    return zone.lights.enduse.power.Mag() * dt;
}

double office::update_plugs(double dt)
{
    return zone.plugs.enduse.power.Mag() * dt;
}

double office::update_hvac(double dt)
{
    const double &Tout = (*(zone.current.pTemperature));
    const double Trange = 40;   /* range over which HP works in heating mode */
    const double Taux = zone.hvac.heating.balance_temperature-Trange;
    const double &Tecon = zone.hvac.cooling.balance_temperature;
    const double &TbalHeat = zone.hvac.heating.balance_temperature;
    const double &TmaxCool = zone.hvac.cooling.design_temperature;
    HCMODE &mode = zone.hvac.mode;

    switch (mode) {
    case HC_OFF:
        cop = 0;
        Qrated = 0;
        break;
    case HC_HEAT:
        cop = 1.0 + (zone.hvac.heating.cop-1)*(Tout-Taux)/Trange;
        Qrated = zone.hvac.heating.capacity + zone.hvac.heating.capacity_perF*(zone.hvac.heating.balance_temperature-Tout);
        break;
    case HC_AUX:
        cop = 1.0;
        Qrated = zone.hvac.heating.capacity;
        break;
    case HC_COOL:
        cop = -1.0 - (zone.hvac.cooling.cop+1)*(Tout-TmaxCool)/(TmaxCool-Tecon);
        Qrated = zone.hvac.cooling.capacity - zone.hvac.cooling.capacity_perF*(Tout-zone.hvac.cooling.balance_temperature);
        break;
    case HC_ECON:
        cop = 0.0;
        Qrated = zone.hvac.cooling.capacity;
        break;
    default:
        throw "hvac mode is invalid";
        break;
    }
    if (cop!=0)
        zone.hvac.enduse.power.SetPowerFactor(Qrated/cop/1000,zone.hvac.enduse.power_factor);
    else
        zone.hvac.enduse.power = complex(0,0);
    if (zone.hvac.enduse.power.Re()<0)
        throw "hvac unit is generating electricity!?!";
    else if (!isfinite(zone.hvac.enduse.power.Re()) || !isfinite(zone.hvac.enduse.power.Im()))
        throw "hvac power is not finite!?!";
    return Qrated*dt;
}

TIMESTAMP office::plc(TIMESTAMP t0, TIMESTAMP t1)
{
#define Tout (*(zone.current.pTemperature))
#define Ti (zone.current.air_temperature)
#define DUTY_CYCLE 0.5 
#define mode (zone.hvac.mode)
#define Taux (zone.hvac.heating.balance_temperature-40)
#define Tecon (zone.hvac.cooling.balance_temperature)
#define DELAY (TS_SECOND*120)
    TIMESTAMP t2=TS_NEVER;

    /* compute the new control temperature */
    update_control_setpoints();

    if (Ti<TheatOn) 
    {
        if (Tout<=Taux)
            mode = HC_AUX;
        else
            mode = HC_HEAT;
    }
    else if (Ti>TheatOff && Ti<TcoolOff) 
        mode = HC_OFF;
    else if (Ti>TcoolOn) 
    {
        if (Tout<=Tecon)
            mode = HC_ECON;
        else
            mode = HC_COOL;
    }

    return TS_NEVER;
}

// IMPLEMENTATION OF CORE LINKAGE

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

EXPORT int init_office(OBJECT *obj, OBJECT *parent) 
{
    try {
        if (obj!=NULL)
            return OBJECTDATA(obj,office)->init(parent);
        return 0;
    } catch (char *msg) {
        gl_error("init_%s(obj=%d;%s): %s", obj->oclass->name, obj->id, obj->name?obj->name:"unnamed", msg);
        return 0;
    }
}

EXPORT TIMESTAMP sync_office(OBJECT *obj, TIMESTAMP t1, PASSCONFIG pass)
{
    try {
        TIMESTAMP t2 = TS_NEVER;
        office *my = OBJECTDATA(obj,office);
        switch (pass) {
        case PC_PRETOPDOWN:
            t2 = my->presync(obj->clock,t1);
            break;
        case PC_BOTTOMUP:
            t2 = my->sync(obj->clock,t1);
            break;
        default:
            throw("invalid pass request");
            break;
        }
        if (pass==clockpass)
            obj->clock = t1;        
        return t2;
    } catch (char *msg) {
        gl_error("sync_office(obj=%d;%s): %s", obj->id, obj->name?obj->name:"unnamed", msg);
        return TS_INVALID; /* halt the clock */
    }
}

EXPORT TIMESTAMP plc_office(OBJECT *obj, TIMESTAMP t1, PASSCONFIG pass)
{
    try {
        if (pass==PC_BOTTOMUP)
            return OBJECTDATA(obj,office)->plc(obj->clock,t1);
        else
            return TS_NEVER;
    } catch (char *msg) {
        gl_error("plc_%s(obj=%d;%s): %s", obj->oclass->name, obj->id, obj->name?obj->name:"unnamed", msg);
        return TS_INVALID;
    }
}

---

GridLAB-D^TM Version 1.0

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