# Lecture 14, Algebraic Modeling Languages

Algebraic Modeling Languages (AML) are high-level computer programming languages for describing and solving high complexity problems for large scale mathematical computation (i.e. large scale optimization type problems).  Their syntax mimics the mathematical notation of optimization problems, which allows one to express optimization problems in a familiar, concise and readable way. 

**AMLs do not directly solve the problem, but they call appropriate external solvers to find the solution.**

Examples of AMLs are
* A Mathematical Programming Language (AMPL),
* General Algebraic Modeling System (GAMS),
* Optimization Programming Language (OPL),
* Advanced Interactive Multidimensional Modeling System (AIMMS), and
* Pyomo.

In addition to the ease of modelling, one of the advantages of AMLs is that you can model the problem once and then solve it with multiple solvers.

## Pyomo

On this course, we use Pyomo as an example of AMLs. Pyomo is a Python-based, open-source optimization modeling language with a diverse set of optimization capabilities.

Pyomo may not be a completely typical AML, because Pyomo's modeling objects are embedded within a full-featured high-level programming language providing a rich set of supporting libraries, which distinguishes Pyomo from other AMLs.

Pyomo supports a wide range of problem types, including:
* Linear programming
* Quadratic programming
* Nonlinear programming
* Mixed-integer linear programming
* Mixed-integer quadratic programming
* Mixed-integer nonlinear programming
* Stochastic programming
* Generalized disjunctive programming
* Differential algebraic equations
* Bilevel programming
* Mathematical programs with equilibrium constraints

# Installing Pyomo

The easiest way to install Pyomo is to call
```
pip install pyomo
```
when pip has been installed on your machine.

## Example 1, linear optimization

Let us start with a very simple linear problem
$$
\begin{align}
\min &\qquad   2x_1+3x_2\\
\text{s.t. }& \qquad 3x_1+4x_2\geq 1\\
& \qquad x_1,x_2\geq 0.
\end{align}
$$


```python
from pyomo.environ import *


model = ConcreteModel()

model.x = Var([1,2], domain=NonNegativeReals) #Non-negative variables x[1] and x[2]

model.OBJ = Objective(expr = 2*model.x[1] + 3*model.x[2]) #Objective function

model.Constraint1 = Constraint(expr = 3*model.x[1] + 4*model.x[2] >= 1) #Constraint

```

Once we have defined the problem, we can solve it. Let us start by using glpk, which is an open source linear programming program.

You need to have glpk installed on your system. For details, see https://www.gnu.org/software/glpk/#TOCdownloading. For many Linux distributions, you can install glpk from the repositories by typing
```
sudo yum install glpk
```
```
sudo apt-get install glpk,
```
or whatever your distribution needs.



```python
from pyomo.opt import SolverFactory #Import interfaces to solvers
opt = SolverFactory("glpk") #Use glpk
res = opt.solve(model, tee=True) #Solve the  problem and print the output
print "Solution:"
print "========="
model.x.display() #Print values of x
```

    GLPSOL: GLPK LP/MIP Solver, v4.52
    Parameter(s) specified in the command line:
     --write /tmp/tmpogS1Ue.glpk.raw --wglp /tmp/tmpA5Hc_F.glpk.glp --cpxlp /tmp/tmpkKjSRP.pyomo.lp
    Reading problem data from '/tmp/tmpkKjSRP.pyomo.lp'...
    2 rows, 3 columns, 3 non-zeros
    21 lines were read
    Writing problem data to `/tmp/tmpA5Hc_F.glpk.glp'...
    15 lines were written
    GLPK Simplex Optimizer, v4.52
    2 rows, 3 columns, 3 non-zeros
    Preprocessing...
    1 row, 2 columns, 2 non-zeros
    Scaling...
     A: min|aij| =  3.000e+00  max|aij| =  4.000e+00  ratio =  1.333e+00
    Problem data seem to be well scaled
    Constructing initial basis...
    Size of triangular part is 1
          0: obj =   0.000000000e+00  infeas =  1.000e+00 (0)
    *     1: obj =   7.500000000e-01  infeas =  0.000e+00 (0)
    *     2: obj =   6.666666667e-01  infeas =  0.000e+00 (0)
    OPTIMAL LP SOLUTION FOUND
    Time used:   0.0 secs
    Memory used: 0.0 Mb (40408 bytes)
    Writing basic solution to `/tmp/tmpogS1Ue.glpk.raw'...
    7 lines were written
    Solution:
    =========
    x : Size=2, Index=x_index, Domain=NonNegativeReals
        Key : Lower : Value          : Upper : Fixed : Stale
          1 :     0 : 0.333333333333 :  None : False : False
          2 :     0 :            0.0 :  None : False : False


Now, if you have other linear solvers installed on your system, you can use them too. Let us use Cplex, which is a commercial solvers (academic license available).


```python
opt = SolverFactory("cplex")
res = opt.solve(model, tee=True)
print "Solution:"
model.x.display()
```

    
    Welcome to IBM(R) ILOG(R) CPLEX(R) Interactive Optimizer 12.6.2.0
      with Simplex, Mixed Integer & Barrier Optimizers
    5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55 5655-Y21
    Copyright IBM Corp. 1988, 2015.  All Rights Reserved.
    
    Type 'help' for a list of available commands.
    Type 'help' followed by a command name for more
    information on commands.
    
    CPLEX> Logfile 'cplex.log' closed.
    Logfile '/tmp/tmpnxw9d9.cplex.log' open.
    CPLEX> Problem '/tmp/tmpEdgLVy.pyomo.lp' read.
    Read time = 0.00 sec. (0.00 ticks)
    CPLEX> Problem name         : /tmp/tmpEdgLVy.pyomo.lp
    Objective sense      : Minimize
    Variables            :       3
    Objective nonzeros   :       2
    Linear constraints   :       2  [Greater: 1,  Equal: 1]
      Nonzeros           :       3
      RHS nonzeros       :       2
    
    Variables            : Min LB: 0.000000         Max UB: all infinite   
    Objective nonzeros   : Min   : 2.000000         Max   : 3.000000       
    Linear constraints   :
      Nonzeros           : Min   : 1.000000         Max   : 4.000000       
      RHS nonzeros       : Min   : 1.000000         Max   : 1.000000       
    CPLEX> Tried aggregator 1 time.
    LP Presolve eliminated 2 rows and 3 columns.
    All rows and columns eliminated.
    Presolve time = 0.00 sec. (0.00 ticks)
    
    Dual simplex - Optimal:  Objective =  6.6666666667e-01
    Solution time =    0.00 sec.  Iterations = 0 (0)
    Deterministic time = 0.00 ticks  (2.72 ticks/sec)
    
    CPLEX> Solution written to file '/tmp/tmpjWjvnD.cplex.sol'.
    CPLEX> Solution:
    x : Size=2, Index=x_index, Domain=NonNegativeReals
        Key : Lower : Value          : Upper : Fixed : Stale
          1 :     0 : 0.333333333333 :  None : False : False
          2 :     0 :            0.0 :  None : False : False


We can use also gurobi, which is another commercial solver with academic license.


```python
opt = SolverFactory("gurobi")
res = opt.solve(model, tee=True)
print "Solution:"
model.x.display()
```

    Optimize a model with 2 rows, 3 columns and 3 nonzeros
    Coefficient statistics:
      Matrix range    [1e+00, 4e+00]
      Objective range [2e+00, 3e+00]
      Bounds range    [0e+00, 0e+00]
      RHS range       [1e+00, 1e+00]
    Presolve removed 2 rows and 3 columns
    Presolve time: 0.00s
    Presolve: All rows and columns removed
    Iteration    Objective       Primal Inf.    Dual Inf.      Time
           0    6.6666667e-01   0.000000e+00   0.000000e+00      0s
    
    Solved in 0 iterations and 0.00 seconds
    Optimal objective  6.666666667e-01
    Solution:
    x : Size=2, Index=x_index, Domain=NonNegativeReals
        Key : Lower : Value          : Upper : Fixed : Stale
          1 :     0 : 0.333333333333 :  None : False : False
          2 :     0 :            0.0 :  None : False : False


## Example 2, nonlinear optimization

Let use define optimization problem
$$
\begin{align}
\max &\qquad c_b\\
\text{s.t. }& \qquad c_{af}s_v - s_vc_a-k_1c_a=0\\
&\qquad s_vc_b+k_1c_a-k_2c_b-2k_3c_a^2=0\\
&\qquad s_vc_c+k_2c_b=0\\
&\qquad s_vc_d+k_3c_a^2=0,\\
&\qquad s_v,c_a,c_b,c_c,c_d\geq0
\end{align}
$$
where $k_1=5/6$, $k_2=5/3$, $k_3=1/6000$, and $c_{af}=10000$.


```python
from pyomo.environ import *
# create the concrete model
model = ConcreteModel()
# set the data 
k1 = 5.0/6.0 
k2 = 5.0/3.0 
k3 = 1.0/6000.0 
caf = 10000.0 
# create the variables
model.sv = Var(initialize = 1.0, within=PositiveReals)
model.ca = Var(initialize = 5000.0, within=PositiveReals)
model.cb = Var(initialize = 2000.0, within=PositiveReals)
model.cc = Var(initialize = 2000.0, within=PositiveReals)
model.cd = Var(initialize = 1000.0, within=PositiveReals)

# create the objective
model.obj = Objective(expr = model.cb, sense=maximize)
# create the constraints
model.ca_bal = Constraint(expr = (0 == model.sv * caf \
    - model.sv * model.ca - k1 * model.ca \
    - 2.0 * k3 * model.ca ** 2.0))
model.cb_bal = Constraint(expr=(0 == -model.sv * model.cb \
    + k1 * model.ca - k2 * model.cb))
model.cc_bal = Constraint(expr=(0 == -model.sv * model.cc \
    + k2 * model.cb))
model.cd_bal = Constraint(expr=(0 == -model.sv * model.cd \
    + k3 * model.ca ** 2.0))
```

## Solving with Ipopt

Install IPopt following http://www.coin-or.org/Ipopt/documentation/node10.html.


```python
opt = SolverFactory("ipopt",solver_io="nl")

opt.solve(model,tee=True)

print "Solution is "
model.sv.display()
model.ca.display()
model.cb.display()
model.cc.display()
model.cd.display()
```

    
    
    ******************************************************************************
    This program contains Ipopt, a library for large-scale nonlinear optimization.
     Ipopt is released as open source code under the Eclipse Public License (EPL).
             For more information visit http://projects.coin-or.org/Ipopt
    ******************************************************************************
    
    This is Ipopt version 3.12, running with linear solver mumps.
    NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
    
    Number of nonzeros in equality constraint Jacobian...:       11
    Number of nonzeros in inequality constraint Jacobian.:        0
    Number of nonzeros in Lagrangian Hessian.............:        5
    
    Total number of variables............................:        5
                         variables with only lower bounds:        5
                    variables with lower and upper bounds:        0
                         variables with only upper bounds:        0
    Total number of equality constraints.................:        4
    Total number of inequality constraints...............:        0
            inequality constraints with only lower bounds:        0
       inequality constraints with lower and upper bounds:        0
            inequality constraints with only upper bounds:        0
    
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
       0 -2.0000000e+03 7.50e+03 6.25e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
       1 -1.0801475e+03 5.54e+02 2.46e+00  -1.0 1.39e+03    -  5.54e-01 1.00e+00h  1
       2 -1.0763574e+03 8.86e+01 1.87e+02  -1.0 3.66e+02    -  9.31e-01 1.00e+00h  1
       3 -1.0727252e+03 1.07e+01 5.26e+00  -1.0 9.35e+01    -  9.71e-01 1.00e+00h  1
       4 -1.0726714e+03 4.01e+00 1.22e-01  -1.0 6.03e+01    -  9.68e-01 1.00e+00h  1
       5 -1.0724371e+03 3.44e-04 2.93e-05  -2.5 4.19e-01    -  1.00e+00 1.00e+00h  1
       6 -1.0724372e+03 1.63e-08 3.93e-09  -3.8 3.85e-03    -  1.00e+00 1.00e+00h  1
       7 -1.0724372e+03 5.68e-11 2.68e-11  -5.7 2.24e-04    -  1.00e+00 1.00e+00h  1
       8 -1.0724372e+03 4.55e-13 2.51e-14  -8.6 2.78e-06    -  1.00e+00 1.00e+00h  1
    
    Number of Iterations....: 8
    
                                       (scaled)                 (unscaled)
    Objective...............:  -1.0724372001086319e+03   -1.0724372001086319e+03
    Dual infeasibility......:   2.5080209371475748e-14    2.5080209371475748e-14
    Constraint violation....:   1.1368683772161604e-14    4.5474735088646412e-13
    Complementarity.........:   2.5059065225787342e-09    2.5059065225787342e-09
    Overall NLP error.......:   2.5059065225787342e-09    2.5059065225787342e-09
    
    
    Number of objective function evaluations             = 9
    Number of objective gradient evaluations             = 9
    Number of equality constraint evaluations            = 9
    Number of inequality constraint evaluations          = 0
    Number of equality constraint Jacobian evaluations   = 9
    Number of inequality constraint Jacobian evaluations = 0
    Number of Lagrangian Hessian evaluations             = 8
    Total CPU secs in IPOPT (w/o function evaluations)   =      0.002
    Total CPU secs in NLP function evaluations           =      0.000
    
    EXIT: Optimal Solution Found.
     
    Ipopt 3.12: Optimal Solution Found
    Solution is 
    sv : Size=1, Index=None, Domain=PositiveReals
        Key  : Lower : Value         : Upper : Fixed : Stale
        None :     0 : 1.34381176107 :  None : False : False
    ca : Size=1, Index=None, Domain=PositiveReals
        Key  : Lower : Value         : Upper : Fixed : Stale
        None :     0 : 3874.25886723 :  None : False : False
    cb : Size=1, Index=None, Domain=PositiveReals
        Key  : Lower : Value         : Upper : Fixed : Stale
        None :     0 : 1072.43720011 :  None : False : False
    cc : Size=1, Index=None, Domain=PositiveReals
        Key  : Lower : Value         : Upper : Fixed : Stale
        None :     0 : 1330.09353341 :  None : False : False
    cd : Size=1, Index=None, Domain=PositiveReals
        Key  : Lower : Value         : Upper : Fixed : Stale
        None :     0 : 1861.60519963 :  None : False : False


# Black-box optimization using scipy.optimize

Often cases, you do not have algebraic formulations of the objective functions, but instead, you have an executable, which gives you the values and you do not know what is happening inside there.

### Example

Executable 'prob4' (which you need to compile before using) includes a script for two variable problem with three inequality constraints (where two of them involve only one variable). The problem is of the form
$$
\min \  f(x)
\\ \text{s.t. }g(x) <= 2
\\          x\geq 0.
$$
The executable reads in a file 'input.txt', which contains variable values on top of each other and outputs a file "output.txt", which contains on top of each other value of f, value of g, value of $x_1$, value of $x_2$, gradient of f and gradient of the contraints.

Let us solve this problem using *scipy.optimize*



```python
import csv
def evaluate_prob4(x):
    with open('input.txt','w') as f:
        f.write('%f\n%f'%(x[0],x[1])) #Write x[0] and x[1] to the input.txt file
    !./prob4 #Execute prob4
    val = []
    with open('output.txt','r') as f: 
        valuereader = csv.reader(f)
        for row in valuereader:
            val.extend([float(i) for i in row])
    f_val = val[0]
    g_val = [0]*3
    g_val[0] = 2-val[1]
    g_val[1]=val[2]
    g_val[2]=val[3]
    grad_f=[val[4],val[5]]
    grad_g = [[0,0],[0,0],[0,0]]
    grad_g[0] = [-val[6],-val[7]]
    grad_g[1] = [val[8],val[9]]
    grad_g[2] = [val[10],val[11]]
    return f_val,g_val,grad_f,grad_g
        
```


```python
import math
evaluate_prob4([2,0.])
```




    (3.0, [-2.0, 2.0, 0.0], [-6.0, 2.0], [[-4.0, -0.0], [1.0, 1.0], [1.0, 1.0]])




```python
from scipy.optimize import minimize

constraint_tuple=(
    {'type':'ineq','fun':lambda x:evaluate_prob4(x)[1][0],\
     'jac':lambda x:evaluate_prob4(x)[3][0]},
    {'type':'ineq','fun':lambda x:evaluate_prob4(x)[1][1],\
     'jac':lambda x:evaluate_prob4(x)[3][1]},
    {'type':'ineq','fun':lambda x:evaluate_prob4(x)[1][2],\
     'jac':lambda x:evaluate_prob4(x)[3][2]}
)
```


```python
minimize(lambda x: evaluate_prob4(x)[0], [0,0], method='SLSQP'
                        , jac=lambda x: evaluate_prob4(x)[2], 
         constraints = constraint_tuple,options = {'disp':True})
```

# Example 3, Nonlinear multiobjective optimization

Let us study optimization problem
$$
\begin{align}
\min \ & \left(\sum_{i=1}^{48}\frac{\sqrt{1+x_i^2}}{v_i},\sum_{i=1}^{48}\left(\left(\frac{x_iv_i}{\sqrt{1+x_i^2}}+v_w\right)^2+\frac{v_i^2}{1+x_i^2}\right)\right., \\
&\qquad\left.\sum_{i=1}^{47}\big|(x_{i+1}-x_i\big|\right)\\
\text{s.t. } & \sum_{i=1}^{j}x_i\leq -1\text{ for all }j=10,11,12,13,14\\
& \left|\sum_{i=1}^{j}x_i\right|\geq 2\text{ for all }j=20,21,22,23,24\\
& \sum_{i=1}^{j}x_i\geq 1\text{ for all }j=30,31,32,33,34\\
&\sum_{i=1}^{48}\frac{\sqrt{1+x_i^2}}{v_i} \leq 5\\
&\sum_{i=1}^{48}x_i=0\\
&-10\leq\sum_{i=1}^{j}x_i\leq10\text{ for all }j=1,\ldots,48
&0\leq v_i\leq 25\text{ for all }i=1,\ldots,48\\
&-10\leq x_i\leq 10\text{ for all }i=1,\ldots,48\\
\end{align}
$$


```python

from pyomo.environ import *
# create the concrete model9
def solve_ach(reference,lb,ub):
    model = ConcreteModel()

    vwind = 5.0
    min_speed = 0.01


    #f1, time used
    def f1(model):
        return sum([sqrt(1+model.y[i]**2)/model.v[i] for i in range(48)])
    #f2, wind drag, directly proportional to square of speed wrt. wind
    def f2(model):
        return sum([((model.y[i]*model.v[i])/sqrt(1+model.y[i]**2)+vwind)**2/
                    +model.v[i]**2*((1+model.y[i])**2) for i in range(48)])
    #f3, maximal course changes
    def f3(model):
        return sum([abs(model.y[i+1]-model.y[i]) for i in range(47)])

    def h1_rule(model,i):
        return sum(model.y[j] for j in range(i))<=-1
    def h2_rule(model,i):
        return abs(sum(model.y[j] for j in range(i)))>=2
    def h3_rule(model,i):
        return sum(model.y[j] for j in range(i))>=1
    def h4_rule(model):
        return sum([sqrt(1+model.y[i]**2)/model.v[i] for i in range(48)])<=25
    def h5_rule(model):
        return sum(model.y[i] for i in range(48))==0

    def f_rule(model):
        return t

    def y_init(model,i):
        if i==0:
            return -1
        if i==18:
            return -1
        if i==24:
            return 1
        if i==25:
            return 1
        if i==26:
            return 1
        if i==34:
            return -1
        return 0
    model.y = Var(range(48),bounds = (-10,10),initialize=y_init)
    model.v = Var(range(48),domain=NonNegativeReals,bounds=(min_speed,25),initialize=25)
    model.t = Var()
    model.h1=Constraint(range(9,14),rule=h1_rule)
    model.h2=Constraint(range(19,24),rule=h2_rule)
    model.h3=Constraint(range(29,34),rule=h3_rule)
    model.h4=Constraint(rule=h4_rule)
    model.h5=Constraint(rule=h5_rule)
    
    def h6_rule(model,i):
        return -10<=sum([model.y[j] for j in range(i)])<=10
    
    model.h6 = Constraint(range(1,48),rule=h6_rule)
    def t_con_f1_rule(model):
        return model.t>=(f1(model)-reference[0]-lb[0])/(ub[0]-lb[0])
    model.t_con_f1 = Constraint(rule = t_con_f1_rule)
    def t_con_f2_rule(model):
        return model.t>=(f2(model)-reference[1]-lb[1])/(ub[1]-lb[1])
    model.t_con_f2 = Constraint(rule = t_con_f2_rule)
    def t_con_f3_rule(model):
        return model.t>=(f3(model)-reference[2]-lb[2])/(ub[2]-lb[2])
    model.t_con_f3 = Constraint(rule = t_con_f3_rule)
    model.f = Objective(expr = model.t+1e-10*(f1(model)+f2(model)+f3(model)))
    tee =False
    opt = SolverFactory("ipopt",solver_io="nl")
    opt.options.max_iter=100000
    #opt.options.constr_viol_tol=0.01
    #opt.options.halt_on_ampl_error = "yes"

    opt.solve(model,tee=tee)
    return [[value(f1(model)),value(f2(model)),value(f3(model))],[model.y,model.v]]

```


```python
lb_ = [0,0,0]
ub_ = [1,1,1]
values =[]
for i in range(3):
    reference = [1e10,1e10,1e10]
    reference[i]=0
    values.append(solve_ach(reference,ub_,lb_)[0])
print values
```

    WARNING - Loading a SolverResults object with a warning status into model=unknown; message from solver=Ipopt 3.12\x3a Maximum Number of Iterations Exceeded.
    WARNING - Loading a SolverResults object with a warning status into model=unknown; message from solver=Ipopt 3.12\x3a Maximum Number of Iterations Exceeded.
    [[25133.93493837363, 759264160.1226324, 860.0], [25.000000247539568, 11188.12151236145, 147.55358533550918], [12673.93579526275, 289156195.6539606, 151.2406436251979]]



```python
lb = [0,0,0]
ub = [1,1,1]
for i in range(3):
    lb[i] = min([values[j][i] for j in range(3)])
    ub[i] = max([values[j][i] for j in range(3)])
print lb
print ub
```


```python
[f,x] = solve_ach([(a+b)/2 for (a,b) in zip(lb,ub)],lb,ub) #Compromise solution
#[f,x] = solve_ach([1e10,1e10,0],lb,ub) #Minimize the third objective

```


```python
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
plt.plot([sum(value(x[0][j]) for j in range(i)) for i in range(49)])
currentAxis = plt.gca()
currentAxis.add_patch(Rectangle((10, -1),4,10))
currentAxis.add_patch(Rectangle((20, -2),4,4))
currentAxis.add_patch(Rectangle((30, -10),4,11))
plt.show()
```

b4d5c746201f9f566f2139c7a9949f6d030a4344
ipynb
Jupyter Notebook
Lecture 10, Algebraic Modeling Languages, especially Pyomo.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2019-04-26T12:46:14.000Z
2021-11-23T03:38:59.000Z
Lecture 10, Algebraic Modeling Languages, especially Pyomo.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
Lecture 10, Algebraic Modeling Languages, especially Pyomo.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2016-01-08T16:28:11.000Z
2021-04-10T05:18:10.000Z
Qwen/Qwen-72B
1. YES
2. YES