residential/solvers.cpp

// solvers.cpp: general equation solver for compound exponential decays
//
// ETP solver utilities
//
// DP Chassin
// December 2002
// Copyright (C) 2002, Pacific Northwest National Laboratory
// All Rights Reserved
//
// Update with advanced solver in June 2008 by D. Chassin

#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <time.h>
#include "gridlabd.h"

#include "solvers.h"

#ifndef WIN32
#define _finite isfinite
#endif
static unsigned long _nan[2] = {0xffffffff,0x7fffffff};
static double NaN = *(double*)&_nan;
#define EVAL(t,a,n,b,m,c) (a*exp(n*t) + b*exp(m*t) + c)

double e2solve( double a,   
                double n,   
                double b,   
                double m,   
                double c,   
                double p,   
                double *e)  
{
    double t=0;
    double f = EVAL(t,a,n,b,m,c);

    // is there an extremum/inflexion to consider
    if (a*b<0)
    {
        // compute the time t and value fi at which it occurs
        double an_bm = -a*n/(b*m);
        double tm = log(an_bm)/(m-n);
        double fm = EVAL(tm,a,n,b,m,c);
        double ti = log(an_bm*n/m)/(m-n);
        double fi = EVAL(ti,a,n,b,m,c);
        if (tm>0) // extremum in domain
        {
            if (f*fm<0) // first solution is in range
                t = 0;
            else if (c*fm<0) // second solution is in range
                t = ti;
            else // no solution is in range
                return NaN;
        }
        else if (tm<0 && ti>0) // no extremum but inflexion in domain
        {
            if (fm*c<0) // solution in range
                t = ti;
            else // no solution in range
                return NaN;
        }
        else if (ti<0) // no extremum or inflexion in domain
        {
            if (fi*c<0) // solution in range
                t = ti;
            else // no solution in range
                return NaN;
        }
        else // no solution possible (includes tm==0 and ti==0)
            return NaN;
    }

    // solve using Newton's method
    int iter = 100;
    if (t!=0) // initial t changed to inflexion point
        double f = EVAL(t,a,n,b,m,c);
    double dfdt = EVAL(t,a*n,n,b*m,m,0); 
    while ( fabs(f)>p && isfinite(t) && iter-->0 )
    {
        t -= f/dfdt;
        f = EVAL(t,a,n,b,m,c);
        dfdt = EVAL(t,a*n,n,b*m,m,0);
    }
    if (iter==0)
    {
        gl_warning("etp::solve(a=%.4f,n=%.4f,b=%.4f,m=%.4f,c=%.4f,prec=%.g) failed to converge",a,n,b,m,c,p);
        return NaN; // failed to catch limit condition above
    }
    if (e!=NULL)
        *e = p/dfdt;
    return t>0?t:NaN;
}

#ifdef OLD_SOLVER
#define SGNP(X,P) ((X)<-(P)?-1:((X)>(P)?+1:0))

int dual_decay_solve(double *ans, double prec, double start, double end, int f, double a, double n, double b, double m, double c)
{
#define EVAL(f,t,a,n,b,m,c) (a*pow(n,f)*exp(n*t) + b*pow(m,f)*exp(m*t) + c)

    // check differential for sign change over range
    double df_start = EVAL(f+1,start,a,n,b,m,c);
    double df_end = EVAL(f+1,end,a,n,b,m,c);

    if(!_finite(df_start) || !_finite(df_end))
        return -1;

    int sa = SGNP(df_start,prec);
    int sb = SGNP(df_end,prec);

    // if sign of first differential changed over range
    if (sa*sb==-1)
    {
        // min/max in start-end range implies multiple solutions
        double t;
        if (!dual_decay_solve(&t, prec, start, end, a*pow(n,f+1),n,b*pow(m,f+1),m,0))

            // error: should have a solution for the inflexion point
            return -1;
        
        // look for solutions right and left of first min/max, keep left result
        int n_right = dual_decay_solve(ans, prec, t, end, f, a,n,b,m,c);
        if (n_right<0)
            return -1; // error; right solution failed
        
        int n_left = dual_decay_solve(ans, prec, start, t, f, a,n,b,m,c);
        if (n_left<0)
            return -1; // error: left solution failed

        return n_left+n_right;
    }
    return dual_decay_solve(ans, prec, start, end, a,n,b,m,c);
#undef EVAL
}

int dual_decay_solve(double *ans, double prec, double start, double end, double a, double n, double b, double m, double c)
{
#define EVAL(t,a,n,b,m,c) (a*exp(n*t) + b*exp(m*t) + c)
    // check for sign change over range
    double f_start = EVAL(start,a,n,b,m,c);
    double f_end = EVAL(end,a,n,b,m,c);

    if (!_finite(f_start) ||!_finite(f_end))
        return -1;

    int sa = SGNP(f_start,prec);
    int sb = SGNP(f_end,prec);

    // if sign hasn't change in range
    if (sa*sb==1)

        // no solution possible
        return 0;

    // check for boundary solutions
    else if (sa==0)
        *ans = start;
    else if (sb==0)
        *ans = end;
    else // opposite signs
    {
        // only one solution exists -- use Newton's method to find it
        double t=(start+end)/2;
        int result = dual_decay_solve(ans,prec,start,t,a,n,b,m,c);
        if (result==0)
            result = dual_decay_solve(ans,prec,t,end,a,n,b,m,c);
        return result;
    }   

    return 1;
}
#endif

#undef SGNP

---

GridLAB-D^TM Version 1.0

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