Set the objective of an optimizer as the standard deviation of the input (Nonlinear optimization using pymo)

Question

I am trying to use pymo for a single objective nonlinear optimization problem.

The objective function is to minimize the variance (or standard deviation) of the input variables following certain constraints (which I was able to do in Excel).

Following is a code example of what I am trying to do

model = pyo.ConcreteModel()

# declare decision variables
model.x1 = pyo.Var(domain=pyo.NonNegativeReals)
model.x2 = pyo.Var(domain=pyo.NonNegativeReals)
model.x3 = pyo.Var(domain=pyo.NonNegativeReals)
model.x4 = pyo.Var(domain=pyo.NonNegativeReals)

# declare objective
from statistics import stdev
model.variance = pyo.Objective(
    expr = stdev([model.x1, model.x2, model.x3, model.x4]),
    sense = pyo.minimize)

# declare constraints
model.max_charging = pyo.Constraint(expr = model.x1 + model.x2 + model.x3 + model.x4 >= 500)
model.max_x1 = pyo.Constraint(expr = model.x1 <= 300)
model.max_x2 = pyo.Constraint(expr = model.x2 <= 200)
model.max_x3 = pyo.Constraint(expr = model.x3 <= 100)
model.max_x4 = pyo.Constraint(expr = model.x4 <= 200)

# solve
pyo.SolverFactory('glpk').solve(model).write()

#print
print("energy_price = ", model.variance())
print(f'Variables = [{model.x1()},{model.x2()},{model.x3()},{model.x4()}]')

The error I get is TypeError: can't convert type 'ScalarVar' to numerator/denominator

The problem seems to be caused by using the stdev function from statistics.

My assumption is that the models variables x1-x4 are yet to have been assigned a value and that is the main issue. However, I am not sure how to approach this?

Answer 1

So I managed to solve this issue and I am including the solution below. But first, there are a couple of points I would like to highlight

1. As @Erwin Kalvelagen mentioned, 'stdev' is nonlinear so a linear solver like 'glpk' would always result in an error. For my problem, 'ipopt' worked fine but be careful as it can perform poorly in some cases.
2. Also, as @Erwin Kalvelagen mentioned, Pyomo does not know about the statistics package. So When you try to use a function from that package (e.g., 'stdev', 'variance', etc.), it will try to evaluate the model variables before the solver assigns them any value and that will cause an error.
3. A pyomo objective function needs an expression. The code sample below shows how to dynamically generate an expression for the variance without using an external package. The code is also agnostic to the number of model variables you have.
4. Using either the variance or the standard deviation will serve the same purpose for my project. I opted for using the variance to avoid calculating its square root (as the standard deviation is the square root of the variance).

Variability Objective Function

import pyomo.environ as pyo
def variabilityRule(model):
    #Calculate mean of all model variables
    for index in model.x:
        mean += model.x[index]
    
    mean = mean/len(model.x)

    #Calculate the difference between each model variable and the mean
    obj_exp = ((model.x[1])-mean) * ((model.x[1])-mean)
    for i in range(2, len(model.x)+1): #note that pyomo variables start from 1 not zero
        obj_exp += ((model.x[i])-mean) * ((model.x[i])-mean)

    #Divide by number of items
    obj_exp = (obj_exp/(len(model.x)))

    return obj_exp

model = pyo.ConcreteModel()
model.objective = pyo.Objective(
    rule = variabilityRule,
    sense = pyo.maximize)

EDIT: Standard Deviation Calculation

You can use the same approach to calculate the standard deviation as I found out. Just multiply the final expression (`obj_exp`) by power 0.5

obj_exp = obj_exp ** 0.5

...

P.S. If you are interested, this is how I dynamically generated my model variables based on an input array

model.x = pyo.VarList(domain=pyo.NonNegativeReals)
for i in range(0, len(input_array)):
    model.x.add()