I want to solve an optimization problem where the objective function is the summation of nonlinear, piecewise functions in the decision variables q_i's such that when a decision variable q_i < 1, the function associated with it is quadratic; otherwise it is a log function. In Pyomo with "ipopt" solver, I tried to define the objective function (as suggested by the answer here) using the "Expr_If" expression. However, when I run the code (see below), the solver indicates that an optimal solution is reached (Message: Ipopt 3.14.9\x3a Optimal Solution Found, Termination condition: optimal). I can print the value of optimal decision variables by running:
for x in model.q: print(model.q[x].value)
but I can not print the optimal value of the objective function when I run:
model.payoff()
such that I get this error "math domain error". This error might suggest that a log function is being evaluated at a negative value, but based on the objective function that I defined this should not happen. Also, I can calculate the value of the objective function by evaluating it (after obtaining the optimal solution) directly using the values of optimal solution q_i's by running the code:
payoff_ = 0
for i in model.P:
if model.q[i].value>=1:
payoff_ += model.Beta[i] * log(model.q[i])
else:
payoff_ += (-0.5)*model.Alpha[i] * (model.q[i]-1)**2
print(payoff_())
Do you know why I am getting the "math domain error" when I run model.payoff()
?
Also, is there any other solver I can use to solve this problem? (not necessarily with python)
My code:
model = ConcreteModel()
#Define the set
model.P = Set(initialize=['P1','P2','P3','P4'])
#Parameters
model.Beta = Param(model.P, initialize = {'P1':1,'P2':1.2,'P3':1.4,'P4':1.6})
model.Alpha = Param(model.P, initialize = {'P1':0.1,'P2':0.2,'P3':0.3,'P4':0.4})
#Variables
model.q = Var(model.P)
#Objective
def Payoff(model):
return sum(Expr_if(IF=model.q[i]>=1, THEN=model.Beta[i] * log(model.q[i]),
ELSE=(-0.5)*model.Alpha[i] * (model.q[i]-1)**2) for i in model.P)
model.payoff = Objective(expr = Payoff, sense = maximize)
#Constraints
def limit(model, i):
return -1.1<= model.q[i]
model.limit = Constraint(model.P, rule = limit)
def balance(model):
return summation(model.q) == 0
model.balance = Constraint(rule = balance)
solver = SolverFactory('ipopt')
solver.solve(model)
#model.pprint()
#model.payoff() <--- printing value of objective function gives an error
#value(model.payoff) <--- printing value of objective function gives an error