Search code examples
linear-programmingpyomo

Pyomo: how to sequence a list in a ConcreteModel


I need to make a complex optimization problem. I have tryed some solution for calculating the average of n highest values but I am not reaching the correct solution. The vector that I want to calculate the average of n highest values is Cenpd_VP_TTL. This average will be addressed to CVaR_TTL that will be used in obj.


Solution

  • As stated in comment, you cannot use comparison or logical operators on the variables during the construction of the model as their values are unknown. You can, however, with a little creativity, usually add some variables (likely binary variables) to work through the logic.

    The example below uses some binary logic to (a) figure out the n-highest values of a variable, (b) compute the average of them by adding another variable. Recall, we want to avoid multiplying variables together because that is a non-linear operation.

    This example hinges on a couple of big-M formulations. If you aren't familiar, you should google it or look for a good explanation in an entry-level LP book.

    You didn't have much context in your example about what was being maximized or minimized and the types of constraints, some modifications to the approach below might be needed. This assumes there is some "downward pressure" on the average, so restricting the contribution of each x value in m.x_contribution by limiting the domain to non-negatives works hand-in-hand with the constraint.

    code:

    import pyomo.environ as pyo
    
    m = pyo.ConcreteModel()
    
    m.I = pyo.Set(initialize=range(5))
    
    num_highs_to_select = 2
    
    # VARS
    m.x = pyo.Var(m.I)  # the values that you want to consider for the avg
    m.high_x = pyo.Var(m.I, domain=pyo.Binary)  # 1 if selected as one of the n highest values
    m.x_contribution = pyo.Var(m.I, domain=pyo.NonNegativeReals)  # the contribution from this variable to the average
    
    
    # OBJ:  minimize the average of the n highest x values...
    high_avg = sum(m.x_contribution[i] for i in m.I) / num_highs_to_select
    m.obj = pyo.Objective(expr=high_avg, sense=pyo.minimize)
    
    # some goofy enforcement of x vals to simulate other model constraints....
    @m.Constraint(m.I)
    def limit_x(m, i):
        return m.x[i] == i**2
    
    # limit the number of selected vals
    m.high_select_limit = pyo.Constraint(expr=pyo.sum_product(m.high_x) == num_highs_to_select)
    
    # now the tricky part, we can use a big-M formulation with a suitably large value of M to pick the highest vals
    M = 100
    @m.Constraint(m.I, m.I)
    def compare(m, i, ii):
        # this basically says that every x[i] must be >= every other x[ii] if x[i] is selected AND the other is not selected
        return m.x[i] >= m.x[ii] - (1 - m.high_x[i] + m.high_x[ii]) * M
    
    # gather the contributions for the x values, if they are selected, again with a big-M constraint
    @m.Constraint(m.I)
    def contribute(m, i):
        return m.x_contribution[i] >= m.x[i] - (1 - m.high_x[i]) * M
    
    solver = pyo.SolverFactory('cbc')
    res = solver.solve(m)
    print(res)
    m.x.display()
    m.high_x.display()
    m.x_contribution.display()
    m.obj.display()
    

    yields:

    Problem: 
    - Name: unknown
      Lower bound: 12.5
      Upper bound: 12.5
      Number of objectives: 1
      Number of constraints: 15
      Number of variables: 9
      Number of binary variables: 5
      Number of integer variables: 5
      Number of nonzeros: 4
      Sense: minimize
    Solver: 
    - Status: ok
      User time: -1.0
      System time: 0.0
      Wallclock time: 0.0
      Termination condition: optimal
      Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
      Statistics: 
        Branch and bound: 
          Number of bounded subproblems: 0
          Number of created subproblems: 0
        Black box: 
          Number of iterations: 4
      Error rc: 0
      Time: 0.008842945098876953
    Solution: 
    - number of solutions: 0
      number of solutions displayed: 0
    
    x : Size=5, Index=I
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :  None :   0.0 :  None : False : False :  Reals
          1 :  None :   1.0 :  None : False : False :  Reals
          2 :  None :   4.0 :  None : False : False :  Reals
          3 :  None :   9.0 :  None : False : False :  Reals
          4 :  None :  16.0 :  None : False : False :  Reals
    high_x : Size=5, Index=I
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :   0.0 :     1 : False : False : Binary
          1 :     0 :   0.0 :     1 : False : False : Binary
          2 :     0 :   0.0 :     1 : False : False : Binary
          3 :     0 :   1.0 :     1 : False : False : Binary
          4 :     0 :   1.0 :     1 : False : False : Binary
    x_contribution : Size=5, Index=I
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :   0.0 :  None : False : False : NonNegativeReals
          1 :     0 :   0.0 :  None : False : False : NonNegativeReals
          2 :     0 :   0.0 :  None : False : False : NonNegativeReals
          3 :     0 :   9.0 :  None : False : False : NonNegativeReals
          4 :     0 :  16.0 :  None : False : False : NonNegativeReals
    obj : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True :  12.5
    [Finished in 276ms]