Currently I am building an economic model in Python GEKKO:

m = GEKKO()

m.time = np.arange(0, 200, 1)

#define the variables

consumption                = m.Var()
tax                        = m.Var()
disposable_income          = m.Var()
(...)


#define the parameters

government_spending                         = m.Param(value=20)
interest_rate                               = m.Param(value=0.05)
tax_rate                                    = m.Param(value=0.2)
(...)


#define lagged values of bond_price and bond_holding

m.delay   (bond_price, bond_price_previous, steps=1)
m.delay   (bond_holding, bond_holding_previous, steps=1)


#define the equations of the model

m.Equation(money_holding             == wealth - bill_holding_household - bond_value)
m.Equation(money_demand              == expected_wealth - bill_holding_household - bond_value)
m.Equation(bill_holding_central_bank == bill_supply - bill_holding_household)

m.Equation(bill_supply.dt()          == (government_spending + coupon_rate * bill_supply + bond_holding) - (tax + coupon_rate * bill_holding_central_bank) - bond_price * bond_holding.dt())
m.Equation(money_supply.dt()         == bill_holding_central_bank.dt()) 
(...)   
m.Equation(coupon_rate               == interest_rate)
m.Equation(bond_return               == 1 / bond_price)    

m.options.IMODE=4
m.solve(disp=True, DIAGLEVEL=2)

The model above is solved successfully by GEKKO. Degrees of freedom (DOF) = 0, and the solutions are as expected.

Now, I want to add three more equations to the behavior of the variable bond_price:

m.Equation(bond_price.dt()           == bond_price_increase_step * decision_increase_bond_price)
m.Equation(decision_increase_bond_price == m.if2(target_proportion - ceiling_target_proportion, 0, 1))
m.Equation(target_proportion         == bond_value / (bond_value + bill_holding_household))

Previously, bond_price is a parameter, but now it is a variable. It will change depending on whether target_proportion is higher than the ceiling_target_proportion, a parameter, or not. If it is, decision_increase_bond_price becomes 1, and makes bond_price increase by bond_price_increase_step.

The issue with these three equations is that by introducing the if2, DOF is now larger than 0. To be precise, DOF is now 398, which is 2 times 199, the number of time periods in the system as defined by m.time. This possibly means that there are two more variables in the system than the number of equations.

This positive DOF is not due to me not defining some variables. also I am extra confident that the DOF is caused by if2, because from my experiment, removing and adding an if2 reduces and increases my DOFs by 2.

I have some questions regarding my situation:

1. Why is if2 adding two extra variables? Does the reason apply to other logical functions in GEKKO (max, min, if3, etc.) ?

2. Can this kind of binary switch be implemented in GEKKO to ensure that DOF = 0?

Here is the full code: (the model should run without the three new equations and their associated variables and parameters)

from gekko import GEKKO

m = GEKKO()

m.time = np.arange(0, 200, 1)

consumption                = m.Var()
tax                        = m.Var()
disposable_income          = m.Var()

money_holding              = m.Var()
money_demand               = m.Var()
money_supply               = m.Var()
income                     = m.Var()
bill_holding_household     = m.Var()
wealth                     = m.Var()

coupon_rate                = m.Var()
bill_supply                = m.Var()
bill_holding_central_bank  = m.Var()
expected_disposable_income = m.Var()
expected_wealth            = m.Var()
bond_holding               = m.Var()
bond_return                = m.Var()
bond_value                 = m.Var()
bond_price_previous        = m.Var()
bond_holding_previous      = m.Var()
capital_gain               = m.Var()

bond_price                                  = m.Var(value=10)
decision_increase_bond_price                = m.Var()
target_proportion                           = m.Var()

government_spending                         = m.Param(value=20)
interest_rate                               = m.Param(value=0.05)
tax_rate                                    = m.Param(value=0.2)
marginal_propensity_to_consume_income       = m.Param(value=0.6)
marginal_propensity_to_consume_wealth       = m.Param(value=0.4)

expected_disposable_income_adjustment_speed = m.Param(value=0.5)
expected_wealth_adjustment_speed            = m.Param(value=0.5)   
portion_household_expect_increased_return   = m.Param(value=0.1)

bond_price_increase_step                    = m.Param(value=0.01)
ceiling_target_proportion                   = m.Param(value=0.505)
floor_target_proportion                     = m.Param(value=0.495)

lambda20                                    = m.Param(value=0.35)
lambda22                                    = m.Param(value=1)
lambda23                                    = m.Param(value=1)
lambda24                                    = m.Param(value=0.03)
lambda30                                    = m.Param(value=0.3)
lambda32                                    = m.Param(value=1)
lambda33                                    = m.Param(value=1)
lambda34                                    = m.Param(value=0.03)

m.delay   (bond_price, bond_price_previous, steps=1)
m.delay   (bond_holding, bond_holding_previous, steps=1)

m.Equation(money_holding             == wealth - bill_holding_household - bond_value)
m.Equation(money_demand              == expected_wealth - bill_holding_household - bond_value)
m.Equation(bill_holding_central_bank == bill_supply - bill_holding_household)

m.Equation(bill_supply.dt()          == (government_spending + coupon_rate * bill_supply + bond_holding) - (tax + coupon_rate * bill_holding_central_bank) - bond_price * bond_holding.dt())
m.Equation(money_supply.dt()         == bill_holding_central_bank.dt())   
m.Equation(wealth.dt()               == disposable_income - consumption + capital_gain)
m.Equation(expected_disposable_income.dt() == expected_disposable_income_adjustment_speed * (disposable_income - expected_disposable_income))
m.Equation(expected_wealth.dt()      == expected_wealth_adjustment_speed * (wealth - expected_wealth) + expected_disposable_income - consumption)
m.Equation(capital_gain              == bond_holding_previous * (bond_price - bond_price_previous))

m.Equation(consumption               == marginal_propensity_to_consume_income * expected_disposable_income + marginal_propensity_to_consume_wealth * wealth)
m.Equation(disposable_income         == income + coupon_rate * bill_holding_household + bond_holding - tax)
m.Equation(income                    == consumption + government_spending)
m.Equation(tax                       == tax_rate * (income + coupon_rate * bill_holding_household + bond_holding))

m.Equation(bill_holding_household    == expected_wealth * (lambda20 + lambda22 * coupon_rate + lambda23 * bond_return + lambda24 * (expected_disposable_income/expected_wealth)))
m.Equation(bond_value                == expected_wealth * (lambda30 + lambda32 * coupon_rate + lambda33 * bond_return + lambda34 * (expected_disposable_income/expected_wealth)))
m.Equation(bond_holding              == bond_value / bond_price)

m.Equation(coupon_rate               == interest_rate)
m.Equation(bond_return               == 1 / bond_price)    

m.Equation(bond_price.dt()           == bond_price_increase_step * (decision_increase_bond_price))
m.Equation(decision_increase_bond_price == m.if2(target_proportion - ceiling_target_proportion, 0, 1))
m.Equation(target_proportion         == bond_value / (bond_value + bill_holding_household))

m.options.IMODE=4
m.solve(disp=True)

I have tried an imperfect solution, which is substituting the if2 with an inequality:decision_increase_bond_price * (ceiling_target_proportion - target_proportion) <= 0. I expected this would remove the DOF issue, and sure it did. But the model returns "Solution Not Found".

While I do not know why the model doesn't solve, I know that with this inequality, decision_increase_bond_price does not need to toggle. It can just stay at 0 to satisfy the inequality.

Currently I am building an economic model in Python GEKKO:

m = GEKKO()

m.time = np.arange(0, 200, 1)

#define the variables

consumption                = m.Var()
tax                        = m.Var()
disposable_income          = m.Var()
(...)


#define the parameters

government_spending                         = m.Param(value=20)
interest_rate                               = m.Param(value=0.05)
tax_rate                                    = m.Param(value=0.2)
(...)


#define lagged values of bond_price and bond_holding

m.delay   (bond_price, bond_price_previous, steps=1)
m.delay   (bond_holding, bond_holding_previous, steps=1)


#define the equations of the model

m.Equation(money_holding             == wealth - bill_holding_household - bond_value)
m.Equation(money_demand              == expected_wealth - bill_holding_household - bond_value)
m.Equation(bill_holding_central_bank == bill_supply - bill_holding_household)

m.Equation(bill_supply.dt()          == (government_spending + coupon_rate * bill_supply + bond_holding) - (tax + coupon_rate * bill_holding_central_bank) - bond_price * bond_holding.dt())
m.Equation(money_supply.dt()         == bill_holding_central_bank.dt()) 
(...)   
m.Equation(coupon_rate               == interest_rate)
m.Equation(bond_return               == 1 / bond_price)    

m.options.IMODE=4
m.solve(disp=True, DIAGLEVEL=2)

The model above is solved successfully by GEKKO. Degrees of freedom (DOF) = 0, and the solutions are as expected.

Now, I want to add three more equations to the behavior of the variable bond_price:

m.Equation(bond_price.dt()           == bond_price_increase_step * decision_increase_bond_price)
m.Equation(decision_increase_bond_price == m.if2(target_proportion - ceiling_target_proportion, 0, 1))
m.Equation(target_proportion         == bond_value / (bond_value + bill_holding_household))

Previously, bond_price is a parameter, but now it is a variable. It will change depending on whether target_proportion is higher than the ceiling_target_proportion, a parameter, or not. If it is, decision_increase_bond_price becomes 1, and makes bond_price increase by bond_price_increase_step.

The issue with these three equations is that by introducing the if2, DOF is now larger than 0. To be precise, DOF is now 398, which is 2 times 199, the number of time periods in the system as defined by m.time. This possibly means that there are two more variables in the system than the number of equations.

This positive DOF is not due to me not defining some variables. also I am extra confident that the DOF is caused by if2, because from my experiment, removing and adding an if2 reduces and increases my DOFs by 2.

I have some questions regarding my situation:

1. Why is if2 adding two extra variables? Does the reason apply to other logical functions in GEKKO (max, min, if3, etc.) ?

2. Can this kind of binary switch be implemented in GEKKO to ensure that DOF = 0?

Here is the full code: (the model should run without the three new equations and their associated variables and parameters)

from gekko import GEKKO

m = GEKKO()

m.time = np.arange(0, 200, 1)

consumption                = m.Var()
tax                        = m.Var()
disposable_income          = m.Var()

money_holding              = m.Var()
money_demand               = m.Var()
money_supply               = m.Var()
income                     = m.Var()
bill_holding_household     = m.Var()
wealth                     = m.Var()

coupon_rate                = m.Var()
bill_supply                = m.Var()
bill_holding_central_bank  = m.Var()
expected_disposable_income = m.Var()
expected_wealth            = m.Var()
bond_holding               = m.Var()
bond_return                = m.Var()
bond_value                 = m.Var()
bond_price_previous        = m.Var()
bond_holding_previous      = m.Var()
capital_gain               = m.Var()

bond_price                                  = m.Var(value=10)
decision_increase_bond_price                = m.Var()
target_proportion                           = m.Var()

government_spending                         = m.Param(value=20)
interest_rate                               = m.Param(value=0.05)
tax_rate                                    = m.Param(value=0.2)
marginal_propensity_to_consume_income       = m.Param(value=0.6)
marginal_propensity_to_consume_wealth       = m.Param(value=0.4)

expected_disposable_income_adjustment_speed = m.Param(value=0.5)
expected_wealth_adjustment_speed            = m.Param(value=0.5)   
portion_household_expect_increased_return   = m.Param(value=0.1)

bond_price_increase_step                    = m.Param(value=0.01)
ceiling_target_proportion                   = m.Param(value=0.505)
floor_target_proportion                     = m.Param(value=0.495)

lambda20                                    = m.Param(value=0.35)
lambda22                                    = m.Param(value=1)
lambda23                                    = m.Param(value=1)
lambda24                                    = m.Param(value=0.03)
lambda30                                    = m.Param(value=0.3)
lambda32                                    = m.Param(value=1)
lambda33                                    = m.Param(value=1)
lambda34                                    = m.Param(value=0.03)

m.delay   (bond_price, bond_price_previous, steps=1)
m.delay   (bond_holding, bond_holding_previous, steps=1)

m.Equation(money_holding             == wealth - bill_holding_household - bond_value)
m.Equation(money_demand              == expected_wealth - bill_holding_household - bond_value)
m.Equation(bill_holding_central_bank == bill_supply - bill_holding_household)

m.Equation(bill_supply.dt()          == (government_spending + coupon_rate * bill_supply + bond_holding) - (tax + coupon_rate * bill_holding_central_bank) - bond_price * bond_holding.dt())
m.Equation(money_supply.dt()         == bill_holding_central_bank.dt())   
m.Equation(wealth.dt()               == disposable_income - consumption + capital_gain)
m.Equation(expected_disposable_income.dt() == expected_disposable_income_adjustment_speed * (disposable_income - expected_disposable_income))
m.Equation(expected_wealth.dt()      == expected_wealth_adjustment_speed * (wealth - expected_wealth) + expected_disposable_income - consumption)
m.Equation(capital_gain              == bond_holding_previous * (bond_price - bond_price_previous))

m.Equation(consumption               == marginal_propensity_to_consume_income * expected_disposable_income + marginal_propensity_to_consume_wealth * wealth)
m.Equation(disposable_income         == income + coupon_rate * bill_holding_household + bond_holding - tax)
m.Equation(income                    == consumption + government_spending)
m.Equation(tax                       == tax_rate * (income + coupon_rate * bill_holding_household + bond_holding))

m.Equation(bill_holding_household    == expected_wealth * (lambda20 + lambda22 * coupon_rate + lambda23 * bond_return + lambda24 * (expected_disposable_income/expected_wealth)))
m.Equation(bond_value                == expected_wealth * (lambda30 + lambda32 * coupon_rate + lambda33 * bond_return + lambda34 * (expected_disposable_income/expected_wealth)))
m.Equation(bond_holding              == bond_value / bond_price)

m.Equation(coupon_rate               == interest_rate)
m.Equation(bond_return               == 1 / bond_price)    

m.Equation(bond_price.dt()           == bond_price_increase_step * (decision_increase_bond_price))
m.Equation(decision_increase_bond_price == m.if2(target_proportion - ceiling_target_proportion, 0, 1))
m.Equation(target_proportion         == bond_value / (bond_value + bill_holding_household))

m.options.IMODE=4
m.solve(disp=True)

I have tried an imperfect solution, which is substituting the if2 with an inequality:decision_increase_bond_price * (ceiling_target_proportion - target_proportion) <= 0. I expected this would remove the DOF issue, and sure it did. But the model returns "Solution Not Found".

While I do not know why the model doesn't solve, I know that with this inequality, decision_increase_bond_price does not need to toggle. It can just stay at 0 to satisfy the inequality.

• python
• if-statement
• binary
• modeling
• gekko

Share Improve this question Follow edited Mar 15 at 4:47 John Hedengren 14.5k11 gold badge2323 silver badges3131 bronze badges asked Mar 13 at 13:00 user29971580user29971580 2122 bronze badges

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 1

Introducing logical conditions like if2 or if3 in gekko requires a switch IMODE to 6 (dynamic optimization) instead of 4 (simulation). This is because logical conditions introduce additional variables and constraints that require optimization techniques to resolve.

The function if2(condition, value_if_condition<0, value_if_condition>0) is an if-else statement that uses a complementarity constraint. It introduces additional internal variables and equations to ensure differentiability. The same applies to other logical functions like abs2, abs3, max2, max3, min2, min3, and if3.

• if2 formed with complementarity constraints and can be solved with a nonlinear programming (NLP) solver.
• if3 uses binary variables and requires a mixed-integer solver and is more reliable but computationally expensive.

To properly implement this binary switch:

1. Switch to IMODE = 6 to allow GEKKO to handle the logic as an optimization problem:

   m.options.IMODE = 6
   

2. If convergence is a problem, consider using if3 instead of if2, because if3 explicitly introduces a binary variable and enforces the switch without a saddle point at the switching condition that is introduced with the if2 complementarity constraint. Here is a complete example of a switching function but with IMODE=2 as a steady-state optimization problem:

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
m = GEKKO(remote=False)
p = m.Param(value=np.linspace(0,10,41))
# conditional statement
y = m.if3(p-4,p**2,-0.2*(p-4)+7)
m.options.IMODE = 2
m.solve()
lbl = r'$y=\mathrm{if3}(p-4,p^2,-0.2(p-4)+7)$'
plt.plot(p,y,'bo',label=lbl)
plt.text(1,5,r'$p^2$')
plt.text(5,10,r'$-0.2 (p-4)+7$')
plt.legend(loc=4)
plt.ylabel('y')
plt.show()

Additional details on complementarity constraints and binary variable forms in the 2 and 3 function forms is available in the Design Optimization course.