# 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 2x_1+3x_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
```

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()
```

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()
```

## Example 2, nonlinear optimization

Let use define optimization problem
$$
\begin{align}
\min &\qquad c_d\\
\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=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()
```

# 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
```

    [25.000000247539568, 11188.12151236145, 147.55358533550918]
    [25133.93493837363, 759264160.1226324, 860.0]



```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()
```

## Black-box optimization problem


```python
import sys
from pyomo.opt.blackbox import RealOptProblem

class RealProblem1(RealOptProblem):

    def __init__(self):
        RealOptProblem.__init__(self)
        self.lower=[0.0, -1.0, 1.0, None]
        self.upper=[None, 0.0, 2.0, -1.0]
        self.nvars=4

    def function_value(self, point):
        self.validate(point)
        return point.vars[0] - point.vars[1] + (point.vars[2]-1.5)**2 + (point.vars[3]+2)**4

problem = RealProblem1()

```

9083d5861c33e8c53bf6d70bbb2937cf5fc14a5d
ipynb
Jupyter Notebook
Lecture 14, Algebraic Modeling Languages, especially Pyomo.ipynb
AwelEshetu/edxCourses
2fb072d8e45d203080d40ca383de14257a1fcbd2
CC-BY-3.0
2016-07-18T12:49:18.000Z
2017-02-21T12:58:16.000Z
Lecture 14, Algebraic Modeling Languages, especially Pyomo.ipynb
AwelEshetu/edxCourses
2fb072d8e45d203080d40ca383de14257a1fcbd2
CC-BY-3.0
Lecture 14, Algebraic Modeling Languages, especially Pyomo.ipynb
AwelEshetu/edxCourses
2fb072d8e45d203080d40ca383de14257a1fcbd2
CC-BY-3.0
Qwen/Qwen-72B
1. YES
2. YES