When solving nonlinear regression problem over a long data set, the Objective value comes high? As a result estimation fails

user_19301709 0
Viewed 35

I am trying to minimize the average root mean squared error of the following form in a class file using Gekko:

objective = np.sqrt((np.sum((ym-np.array(y))**2))/N/
                    (np.sum((ym-np.mean(ym))**2))/N)

Here is the code:

# Code 
from math import ceil
import numpy as np
import pandas as pd
import os
from gekko import GEKKO
import sys
from demandlib.tools import add_weekdays2df
import matplotlib.pyplot as plt

class HeatBuilding_Personalized:
    
    
    def __init__(self, df_index, **kwargs):
        self.datapath = kwargs.get(
            'datapath', os.path.join(os.path.dirname(__file__), 'bdew_data'))
        self.df = pd.DataFrame(index=df_index)
        self.df = add_weekdays2df(self.df, holiday_is_sunday=True,
                                  holidays=kwargs.get('holidays'))
        self.df['hour'] = self.df.index.hour + 1  # hour of the day
        self.temperature = kwargs.get('temperature')
        self.annual_heat_demand = kwargs.get('annual_heat_demand')
        self.shlp_type = kwargs.get('shlp_type').upper()
        self.wind_class = kwargs.get('wind_class')
        self.building_class = kwargs.get('building_class', 0)
        self.ww_incl = kwargs.get('ww_incl', True)
        self.name = kwargs.get('name', self.shlp_type)
        self.data_points = kwargs.get('data_points')
        self.st_p= kwargs.get('st_p')
        self.end= kwargs.get('end')

    def get_bdew_profile(self):

        """ Calculation of the normalized hourly heat demand
        """
        self.df['temperature'] = self.temperature.values
        self.df['temperature_geo'] = self.weighted_temperature(how='geometric_series')

        sf = self.get_sf_values()

          
        f = self.get_weekday_parameters()
        
        # measurements
        self.df['data_points']=self.data_points.values
        self.df= self.df[self.st_p:self.end]
        self.df=self.df.dropna()
        self.annual_heat_demand= self.df['data_points'].sum()
        self.temperature= pd.DataFrame(self.df['temperature'])
        print(self.df)
        ym = pd.DataFrame(self.df['data_points'])
        print("amount of nan",str(ym.isnull().sum()))
        ymeas_mean = np.mean(ym)
        print(ym)
        print('ymeas_mean:',ymeas_mean)
        
        x1= np.array(self.df['temperature_geo'])
        x2= np.array(self.get_weekday_parameters())
        x3= np.int(self.annual_heat_demand)
        x4= np.array(self.get_sf_values())
        ym= np.array(ym)
        # GEKKO model
        m = GEKKO(remote=False)       
   
        a = m.FV( 3.7,lb=1,ub=4)
        a.STATUS=1
        b = m.FV(-35.1,lb=-40,ub=-30)
        b.STATUS=1
        c = m.FV(7.1,lb=5,ub=9)
        c.STATUS=1
        d = m.FV( 0.9,lb=0.1,ub=1.5)
        d.STATUS=1
   
        
        # variables
        T_g= m.Param(value=x1)
        f=m.Param(value=x2)
        annual_demand=m.Param(value=x3)
        sf=m.Param(value=x4)
        ymeas = m.Param(value=ym)
        N = len(ym)
        print('index n:',N)
        
        yest = m.CV(value=0)
        yest.FSTATUS=1
        # y=m.Var() # I am defining my state variabel

        # y = m.Var()
        # z.FSTATUS=1
        
        # regression equation
        k = m.Intermediate((a / (1 + (b / (T_g - 40)) ** c) + d))
       
        s=m.Intermediate(np.sum(k*f))
        
        kw=m.Intermediate( 1.0 / (s / 24))
        m.Equation(yest == (k* kw  * f * sf) * annual_demand)
        # objectives
      
        # m.Minimize(((yest-ymeas)/ymeas)**2)
        
        m.Obj(m.sqrt((np.sum((ymeas-yest)**2))/N/(np.sum((ymeas-np.mean(ymeas))**2))/N)) 
        # print('Obj init value = ' + str(object_af.value))
        
        # regression mode
        m.options.IMODE = 2
        m.options.SOLVER = 1 # considering APOPT solver for 1 and IPOPT for 3
       
       
        # optimize
        m.options.MAX_ITER = 20
        m.options.OTOL = 1.0e-10
        m.options.RTOL = 1.0e-10
        m.solve(disp=True)
        # print parameters
        # print('Obj after value = ' + str(vd.value))
        print('Optimized, a = ' + str(a.value[0]))
        print('Optimized, b = ' + str(b.value[0]))
        print('Optimized, c = ' + str(c.value[0]))
        print('Optimized, d = ' + str(d.value[0]))
        # print('Optimized, h = ' + str(h.value))
        # sys.exit()

        print("optimization is ok")
       
        sf = self.get_sf_values()   
        f = self.get_weekday_parameters()
        
        h = (a.value[0] / (1 + (b.value[0] / (self.df['temperature_geo'] - 40)) ** c.value[0]) + d.value[0])
        kw = 1.0 / (sum(h * f) / 24)    #1.0 instead of annual heat demand because the #annual heat demand is already multiplied in get_bdew_profile and divide by 24 to get #daily value
        y = (kw * h * f * sf) * self.annual_heat_demand
        
        objective= np.sqrt((np.sum((ym-np.array(y))**2))/N/(np.sum((ym-np.mean(ym))**2))/N)
        print('objective calculated without Gekko:',objective)
        
        return y

It returns this output:

Outputs: Solver : 
APOPT (v1.0) Solution time : 27.2771999999968 sec 
Objective : 40884011.5968099
Successful solution --------------------------------------------------
Optimized, a = 3.8708321781 Optimized, b = -31.844822393 
Optimized, c = 7.8648564579 Optimized, d = 1.0244814518

The objective value is high. Without Gekko the objective is calculated as 0.01904060781034217. Why is it different?

gekko python
2 Answers

It is hard to diagnose the problem because it is missing the data to run and verify. Here are a couple things to change:

  1. Set yest as a Variable instead of a CV. It is automatically adding squared error terms when you declare a CV with FSTATUS=1 as described in the Dynamic Optimization course and in the Estimator Tuning Lesson. Because you are defining a custom objective, there is no need to declare a CV.
        #yest = m.CV(value=0)
        #yest.FSTATUS=1
        yest = m.Var(value=0)
  1. For a direct comparison, try declaring Intermediate variables to inspect the parts of the objective. Also, use m.sum() instead of np.sum() for the Gekko version of summation. The den piece is a number so it can be pre-calculated before the objective function is defined with ym.
den = (np.sum((ym-np.mean(ym))**2)
m.Obj(m.sqrt((np.sum((ymeas-yest)**2))/(den*N*N)))

Please post complete minimal, verifiable code for more specific help.

user_23669410

so we corrected the code like said in the previous comment but we get problems when we try to take a sum of a gekko parameter and we don't understand why. s is the variable where we calculate a sum and try to use it in the next equation but it doesn't work. Even with numpy.sum() it doesn't sum it up. We get this error: TypeError: x must be a python list of GEKKO parameters, variables, or expressions

Any idea what we need to change on factor k or f so the sum can be achieved?

    x1= np.array(self.df['temperature_geo'])
    x2= np.array(self.get_weekday_parameters())
    x3= np.int(self.annual_heat_demand)
    x4= np.array(self.get_sf_values())
    ym= np.array(ym)
    # GEKKO model
    m = GEKKO(remote=False)       
    
    # variables
    T_g= m.Param(value=x1)
    f=m.Param(value=x2)
    annual_demand=m.Param(value=x3)
    sf=m.Param(value=x4)
    ymeas = m.Param(value=ym)
    
    yest = m.Var(value=0)

    # regression equation
    k = m.Intermediate((a / (1 + (b / (T_g - 40)) ** c) + d))
   
    s=m.Intermediate(m.sum(k*f))
    kw=m.Intermediate( 1.0 / (s / 24))
    m.Equation(yest == (k* kw  * f * sf) * annual_demand)
    
    den = (m.sum((ymeas-np.mean(ymeas))**2))
    ben = m.sum((ymeas-yest)**2)
    m.Obj(m.sqrt((ben)/N/(den)/N)) 
    # print('Obj init value = ' + str(object_af.value))
    
    # regression mode
    m.options.IMODE = 2
    m.options.SOLVER = 1 # considering APOPT solver for 1 and IPOPT for 3
   
   
    # optimize
    m.options.MAX_ITER = 20
    m.options.OTOL = 1.0e-10
    m.options.RTOL = 1.0e-10
    m.solve(disp=True)
user_197153300
Related