Search code examples
python-3.xoptimizationgekko

Negative degrees of fredom when using GEKKO python

enter image description here

I'm trying to solve the optimization problem as above. And my code is as belows.

It worked, but I got the negative degrees of freedom problem.

And the objective value was also negative, which I did not expect to be. I expected the positive one.

I can't understand why this happened and don't know how this problem can be solved.

Can somebody give me a suggestion?

Code

# Import package
from gekko import GEKKO
import numpy as np

# Define parameters
P_CO = 600                      # $/tonCO
beta_CO2 = 1                    # no unit
P_CO2 = 60                      # $/tonCO2eq
E_ref = 3.1022616               # tonCO2eq/tonCO
E_dir = -1.600570692            # tonCO2eq/tonCO
E_indir_others = 0.3339226804   # tonCO2eq/tonCO
E_indir_elec_cons = 18.46607256 # GJ/tonCO
C1_CAPEX = 285695               # no unit
C2_CAPEX = 188.42               # no unit
C1_FOX = 82282                  # no unit
C2_FOX = 24.094                 # no unit
C1_ROX = 4471.5                 # no unit
C2_ROX = 96.034                 # no unit
C1_UOX = 1983.7                 # no unit
C2_UOX = 249.79                 # no unit
r = 0.08                        # discount rate
N = 10                          # number of scenarios
T = 30                          # total time period
GWP_init = 0.338723235          # 2020 Electricity GWP in EU 27 countries
theta_max = 1600000             # Max capacity

# Function to make GWP_EU matrix (TxN matrix)
def Electricity_GWP(GWP_init, n_years, num_episodes):

    GWP_mean = 0.36258224*np.exp(-0.16395611*np.arange(1, n_years+2)) + 0.03091272
    GWP_mean = GWP_mean.reshape(-1,1)
    GWP_Yearly = np.tile(GWP_mean, num_episodes) 

    noise = np.zeros((n_years+1, num_episodes))
    stdev2050 = GWP_mean[-1] * 0.25 
    stdev = np.arange(0, stdev2050 * (1 + 1/n_years), stdev2050/n_years)

    for i in range(n_years+1):
        noise[i,:] = np.random.normal(0, stdev[i], num_episodes) 

    GWP_forecast = GWP_Yearly + noise 

    return GWP_forecast

GWP_EU = Electricity_GWP(GWP_init, T, N) # (T+1)*N matrix
GWP_EU = GWP_EU[1:,:] # T*N matrix

print(np.shape(GWP_EU))

# Build Gekko model
m = GEKKO(remote=False)
theta = m.Array(m.Var, N, lb=0, ub=theta_max)
demand = np.ones((T,1))
demand[0] = 8031887.589
for k in range(1,11):
    demand[k] = demand[k-1] * 1.026 
for k in range(11,21):
    demand[k] = demand[k-1] * 1.016
for k in range(21,T):
    demand[k] = demand[k-1] * 1.011 
demand = 0.12 * demand
demand = np.tile(demand, N) # T*N matrix

print(np.shape(demand))

obj = m.sum([m.sum([((1/(1+r))**(t+1))*((P_CO*m.min3(demand[t,s], theta[s])) \
            + (beta_CO2*P_CO2*m.min3(demand[t,s], theta[s])*(E_ref-E_dir-E_indir_others-E_indir_elec_cons*GWP_EU[t,s])) \
            - (C1_CAPEX+C2_CAPEX*theta[s]+C1_FOX+C2_FOX*theta[s])-(C1_ROX+C2_ROX*m.min3(demand[t,s], theta[s])+C1_UOX+C2_UOX*m.min3(demand[t,s], theta[s]))) for t in range(T)]) for s in range(N)])
m.Maximize(obj/N)
m.solve()

Output message

(30, 10)
(30, 10)
 ----------------------------------------------------------------
 APMonitor, Version 1.0.0
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :  11
   Constants    :  0
   Variables    :  5121
   Intermediates:  0
   Connections  :  321
   Equations    :  3901
   Residuals    :  3901
 
 Number of state variables:    5121
 Number of total equations: -  3911
 Number of slack variables: -  2400
 ---------------------------------------
 Degrees of freedom       :    -1190
 
 * Warning: DOF <= 0
 ----------------------------------------------
 Steady State Optimization with APOPT Solver
 ----------------------------------------------
Iter:     1 I:  0 Tm:     18.61 NLPi:    5 Dpth:    0 Lvs:    0 Obj: -1.87E+09 Gap:  0.00E+00
 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  18.619200000000003 sec
 Objective      :  -1.8677021320161405E+9
 Successful solution
 ---------------------------------------------------
Solution

  • The negative DOF warning is because of the slack variables that are created when using the min3() function. It is only a warning that if all of the inequalities are active then this could lead to an over-specified system of equations (more equations than variables). If there is a successful solution then this warning can be ignored.

    The negative objective is because most solvers require a minimization of the objective. Gekko automatically converts m.Maximize(obj) to m.Minimize(-obj). This is an equivalent objective. If you'd like to report the maximization and the positive objective, use the following at the end:

    print('Objective: ',-m.options.OBJFCNVAL)