Differences
This shows you the differences between two versions of the page.
Both sides previous revision Previous revision Next revision | Previous revision | ||
personal:blog:2017:0203_jump_for_gams_users [2017/02/03 14:45] antonello |
personal:blog:2017:0203_jump_for_gams_users [2023/12/22 11:39] (current) antonello [Further help] |
||
---|---|---|---|
Line 15: | Line 15: | ||
You have plenty of development environment to choose from (e.g. Jupiter, Juno), a clear modern language, the possibility to interface your model with third party libraries.. all of this basically for free.\\ | You have plenty of development environment to choose from (e.g. Jupiter, Juno), a clear modern language, the possibility to interface your model with third party libraries.. all of this basically for free.\\ | ||
It is also, at least for my user case, much faster than GAMS. Aside the preparation of the model to pass to the solver, where it is roughly equivalent, in the solver execution I can benefit of having on my system a version of IPOPT compiled with the much more performing ma27 linear solver, while for GAMS I would have to rely on the embedded version that is compiled with the MUMPS linear solver. That's part of the flexibility you gain in using JuMP in place of GAMS. | It is also, at least for my user case, much faster than GAMS. Aside the preparation of the model to pass to the solver, where it is roughly equivalent, in the solver execution I can benefit of having on my system a version of IPOPT compiled with the much more performing ma27 linear solver, while for GAMS I would have to rely on the embedded version that is compiled with the MUMPS linear solver. That's part of the flexibility you gain in using JuMP in place of GAMS. | ||
- | That's said, for people that don't need such flexibility, | + | That's said, for people that don't need such flexibility, |
Line 31: | Line 31: | ||
Run, only once, the following code to install JuMP language and a couple of open source solvers: | Run, only once, the following code to install JuMP language and a couple of open source solvers: | ||
<code julia> | <code julia> | ||
- | Pkg.update() | + | using Pkg # Load the package manager |
- | Pkg.add(" | + | Pkg.update() |
- | Pkg.add(" | + | Pkg.add(" |
- | Pkg.add(" | + | Pkg.add(" |
- | Pkg.add(" | + | Pkg.add(" |
+ | Pkg.add(" | ||
+ | Pkg.add(" | ||
</ | </ | ||
Line 42: | Line 44: | ||
==== Importing the libraries ==== | ==== Importing the libraries ==== | ||
- | You will need to import as a minima the '' | + | You will need to import as a minima the '' |
- | <code julia> | + | < |
- | # Import of the JuMP and DataFrames modules (the latter | + | # Import of the JuMP, GLPK, CSV and DataFrames modules (the latter |
- | using JuMP, DataFrames | + | using CSV, DataFrames, GLPK, JuMP |
</ | </ | ||
Line 56: | Line 58: | ||
<code julia> | <code julia> | ||
- | # Sets | + | # Define sets # |
+ | # | ||
+ | # | ||
+ | # | ||
plants | plants | ||
markets = [" | markets = [" | ||
Line 64: | Line 69: | ||
==== Definition of the " | ==== Definition of the " | ||
- | Capacity of plants and demand of markets are directly defined as dictionaries, | + | Capacity of plants and demand of markets are directly defined as dictionaries, |
<code julia> | <code julia> | ||
- | # Parameters | + | # Define parameters # |
+ | # Parameters | ||
+ | # | ||
+ | # / | ||
+ | # san-diego | ||
a = Dict( # capacity of plant i in cases | a = Dict( # capacity of plant i in cases | ||
" | " | ||
" | " | ||
) | ) | ||
+ | |||
+ | # | ||
+ | # / | ||
+ | # chicago | ||
+ | # topeka | ||
b = Dict( # demand at market j in cases | b = Dict( # demand at market j in cases | ||
" | " | ||
Line 78: | Line 92: | ||
) | ) | ||
- | # distance in thousands of miles | + | # Table d(i, |
- | d_table = wsv""" | + | # new-york |
+ | # seattle | ||
+ | # san-diego | ||
+ | d_table = CSV.read(IOBuffer(""" | ||
plants | plants | ||
seattle | seattle | ||
san_diego | san_diego | ||
- | """ | + | """ |
d = Dict( (r[: | d = Dict( (r[: | ||
+ | # Here we are converting the table in a " | ||
# r[: | # r[: | ||
# m: the second key | # m: the second key | ||
# r[Symbol(m)]: | # r[Symbol(m)]: | ||
+ | # Scalar f freight in dollars per case per thousand miles /90/ ; | ||
f = 90 # freight in dollars per case per thousand miles | f = 90 # freight in dollars per case per thousand miles | ||
+ | # Parameter c(i, | ||
+ | # c(i,j) = f * d(i,j) / 1000 ; | ||
# We first declare an empty dictionary and then we fill it with the values | # We first declare an empty dictionary and then we fill it with the values | ||
c = Dict() # transport cost in thousands of dollars per case ; | c = Dict() # transport cost in thousands of dollars per case ; | ||
[ c[p,m] = f * d[p,m] / 1000 for p in plants, m in markets] | [ c[p,m] = f * d[p,m] / 1000 for p in plants, m in markets] | ||
</ | </ | ||
+ | The above code take advantage of [[http:// | ||
+ | If we take the creation of the d dictionary as example, without List Comprehensions we would have had to write a nested for loop like: | ||
+ | <code julia> | ||
+ | d = Dict() | ||
+ | for r in eachrow(d_table) | ||
+ | for m in markets | ||
+ | d = (r[: | ||
+ | end | ||
+ | end | ||
+ | </ | ||
+ | Using List Comprehension is however quicker to code and more readable. | ||
Line 100: | Line 132: | ||
Here we declare a JuML optimisation model and we give it a name. This name will be then passed as first argument to all the subsequent operations, like creation of variables, constraints and objective function.\\ | Here we declare a JuML optimisation model and we give it a name. This name will be then passed as first argument to all the subsequent operations, like creation of variables, constraints and objective function.\\ | ||
- | We can, if we wish, works with several models at the same time.\\ | + | The solver engine to use is given as argument of the '' |
- | If we do not specify a solver, we let JuML use a suitable solver for the type of problem. Aside to specify the solver, we can also pass it solver-level options, e.g.: | + | We could pass solver-specific |
- | '' | + | '' |
<code julia> | <code julia> | ||
- | # Model declaration | + | # Model declaration |
- | trmodel = Model() | + | trmodel = Model(GLPK.Optimizer) |
</ | </ | ||
==== Declaration of the model variables ==== | ==== Declaration of the model variables ==== | ||
- | Variables can have multiple-dimensions - that is, being indexed under several indexes -, and bounds are given at the same time as their declaration. | + | Variables can have multiple-dimensions - that is, being indexed under several indexes -, and bounds are given at the same time as their declaration.\\ |
+ | Differently from GAMS, we don't need to define the variable that is on the left hand side of the objective function. | ||
<code julia> | <code julia> | ||
- | # Variables | + | ## Define variables ## |
+ | # | ||
+ | # | ||
+ | # | ||
+ | # Positive Variable x ; | ||
@variables trmodel begin | @variables trmodel begin | ||
x[p in plants, m in markets] >= 0 # shipment quantities in cases | x[p in plants, m in markets] >= 0 # shipment quantities in cases | ||
Line 125: | Line 162: | ||
<code julia> | <code julia> | ||
- | # Constraints | + | ## Define contrains ## |
+ | # supply(i) | ||
+ | # supply(i) .. sum (j, x(i,j)) =l= a(i) | ||
+ | # demand(j) | ||
+ | # demand(j) .. sum(i, x(i,j)) =g= b(j); | ||
@constraints trmodel begin | @constraints trmodel begin | ||
supply[p in plants], | supply[p in plants], | ||
Line 151: | Line 192: | ||
<code julia> | <code julia> | ||
print(trmodel) | print(trmodel) | ||
- | </code | + | </code> |
==== Resolution of the model ==== | ==== Resolution of the model ==== | ||
- | It is at this point that the solver is called and the model is passed to the solver engine for its solution. The return value is the status of the optimisation (": | + | It is at this point that the solver is called and the model is passed to the solver engine for its solution. The return value is the status of the optimisation ('' |
<code julia> | <code julia> | ||
- | status = solve(trmodel) | + | optimize!(trmodel) |
- | </julia> | + | status = termination_status(trmodel) |
+ | </code> | ||
==== Visualisation of the results ==== | ==== Visualisation of the results ==== | ||
+ | While you can do any fancy output you may wish after you retrieve the optimal value of the variables with '' | ||
+ | Notice that you can also easily retrieve the dual value associated to the constraint with '' | ||
<code julia> | <code julia> | ||
- | if status == :Optimal | + | if status == MOI.OPTIMAL |
- | println(" | + | println(" |
- | println(getvalue(x)) | + | println(" |
+ | println(value.(x)) | ||
+ | println(" | ||
+ | [println(" | ||
+ | println(" | ||
+ | [println(" | ||
+ | |||
else | else | ||
println(" | println(" | ||
Line 173: | Line 223: | ||
</ | </ | ||
- | ==== Set definition ==== | ||
- | Sets are created as attributes object of the main model objects and all the information is given as parameter in the constructor function. Specifically, | ||
- | <code python> | ||
- | ## Define sets ## | ||
- | # Sets | ||
- | # | ||
- | # | ||
- | model.i = Set(initialize=[' | ||
- | model.j = Set(initialize=[' | ||
- | </ | ||
- | |||
- | ==== Parameters ==== | ||
- | Parameter objects are created specifying the sets over which they are defined and are initialised with either a python dictionary or a scalar: | ||
- | <code python> | ||
- | ## Define parameters ## | ||
- | # | ||
- | # | ||
- | # / | ||
- | # san-diego | ||
- | # | ||
- | # / | ||
- | # chicago | ||
- | # topeka | ||
- | model.a = Param(model.i, | ||
- | model.b = Param(model.j, | ||
- | # Table d(i, | ||
- | # new-york | ||
- | # seattle | ||
- | # san-diego | ||
- | dtab = { | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | (' | ||
- | } | ||
- | model.d = Param(model.i, | ||
- | # Scalar f freight in dollars per case per thousand miles /90/ ; | ||
- | model.f = Param(initialize=90, | ||
- | </ | ||
- | A third, powerful way to initialize a parameter is using a user-defined function.\\ | ||
- | This function will be automatically called by pyomo with any possible (i,j) set. In this case pyomo will actually call c_init() six times in order to initialize the model.c parameter. | ||
- | <code python> | ||
- | # Parameter c(i, | ||
- | # c(i,j) = f * d(i,j) / 1000 ; | ||
- | def c_init(model, | ||
- | return model.f * model.d[i, | ||
- | model.c = Param(model.i, | ||
- | </ | ||
- | |||
- | ==== Variables ==== | ||
- | Similar to parameters, variables are created specifying their domain(s). For variables we can also specify the upper/lower bounds in the constructor.\\ | ||
- | Differently from GAMS, we don't need to define the variable that is on the left hand side of the objective function. | ||
- | <code python> | ||
- | ## Define variables ## | ||
- | # Variables | ||
- | # | ||
- | # | ||
- | # Positive Variable x ; | ||
- | model.x = Var(model.i, | ||
- | </ | ||
- | |||
- | ==== Constrains ==== | ||
- | At this point, it should not be a surprise that constrains are again defined as model objects with the required information passed as parameter in the constructor function. | ||
- | <code python> | ||
- | ## Define contrains ## | ||
- | # supply(i) | ||
- | # supply(i) .. sum (j, x(i,j)) =l= a(i) | ||
- | def supply_rule(model, | ||
- | return sum(model.x[i, | ||
- | model.supply = Constraint(model.i, | ||
- | # demand(j) | ||
- | # demand(j) .. sum(i, x(i,j)) =g= b(j); | ||
- | def demand_rule(model, | ||
- | return sum(model.x[i, | ||
- | model.demand = Constraint(model.j, | ||
- | </ | ||
- | The above code take advantage of [[https:// | ||
- | If we take the supply_rule as example, this is actually called two times by pyomo (once for each of the elements of i). Without List Comprehensions we would have had to write our function using a for loop, like: | ||
- | <code python> | ||
- | def supply_rule(model, | ||
- | supply = 0.0 | ||
- | for j in model.j: | ||
- | supply += model.x[i, | ||
- | return supply <= model.a[i] | ||
- | </ | ||
- | Using List Comprehension is however quicker to code and more readable. | ||
- | |||
- | ==== Objective & solving ==== | ||
- | The definition of the objective is similar to those of the constrains, except that most solvers require a scalar objective function, hence a unique function, and we can specify the sense (direction) of the optimisation. | ||
- | <code python> | ||
- | ## Define Objective and solve ## | ||
- | # cost define objective function | ||
- | # cost .. z =e= sum((i,j), c(i, | ||
- | # Model transport /all/ ; | ||
- | # Solve transport using lp minimizing z ; | ||
- | def objective_rule(model): | ||
- | return sum(model.c[i, | ||
- | model.objective = Objective(rule=objective_rule, | ||
- | </ | ||
- | As we are here looping over two distinct sets, we can see how List Comprehension really simplifies the code. The objective function could have being written without List Comprehension as: | ||
- | <code python> | ||
- | def objective_rule(model): | ||
- | obj = 0.0 | ||
- | for ki in model.i: | ||
- | for kj in model.j: | ||
- | obj += model.c[ki, | ||
- | return obj | ||
- | </ | ||
- | |||
- | ==== Retrieving the output ==== | ||
- | To retrieve the output and do something with it (either to just display it -like we do here-, to plot a graph with [[http:// | ||
- | This function is called by pyomo after the solver has finished. | ||
- | <code python> | ||
- | ## Display of the output ## | ||
- | # Display x.l, x.m ; | ||
- | def pyomo_postprocess(options=None, | ||
- | model.x.display() | ||
- | </ | ||
- | We can print model structure information with '' | ||
- | Results are also by default saved in a results.json file or, if PyYAML is installed in the system, in results.yml. | ||
==== Editing and running the script ==== | ==== Editing and running the script ==== | ||
- | Differently from GAMS you can use whatever editor environment you wish to code a pyomo script. If you don't need debugging features, a simple text editor like Notepad++ (in windows), gedit or kate (in Linux) will suffice. They already have syntax highlight for python.\\ | + | Differently from GAMS you can use whatever editor environment you wish to code a JuMP script. If you don't need debugging features, a simple text editor like Notepad++ (in windows), gedit or kate (in Linux) will suffice. They already have syntax highlight for Julia.\\ |
- | If you want advanced features and debugging capabilities you can use a dedicated | + | If you want advanced features and debugging capabilities you can use a dedicated |
- | You will normally run the script as '' | + | If you are using instead |
- | If you want to run the script as '' | + | |
- | <code python> | + | |
- | # This is an optional code path that allows the script to be run outside of | + | |
- | # pyomo command-line. For example: | + | |
- | if __name__ == ' | + | |
- | #This replicates what the pyomo command-line tools does | + | |
- | from pyomo.opt import SolverFactory | + | |
- | import pyomo.environ | + | |
- | opt = SolverFactory(" | + | |
- | instance = model.create() | + | |
- | results = opt.solve(instance) | + | |
- | #sends results to stdout | + | |
- | results.write() | + | |
- | pyomo_postprocess(None, | + | |
- | </ | + | |
- | + | ||
- | Finally, if you are very lazy and want to run the script with just ./ | + | |
- | <code python> | + | |
- | # | + | |
- | # -*- coding: utf-8 -*- | + | |
- | </ | + | |
===== Further help ===== | ===== Further help ===== | ||
- | Documentation of pyomo is available from [[https://software.sandia.gov/ | + | Documentation of JuMP is available from [[https://jump.dev/|this page]], |
+ | |||
+ | Happy modelling with JuMP ;-) | ||
- | Happy modelling with pyomo ;-) | ||
===== Complete script ===== | ===== Complete script ===== | ||
Here is the complete script: | Here is the complete script: | ||
- | < | + | < |
- | #!/ | + | # Transport example |
- | # -*- coding: utf-8 -*- | + | |
+ | # Transposition in JuMP of the basic transport model used in the GAMS tutorial | ||
+ | # | ||
+ | # This problem finds a least cost shipping schedule that meets | ||
+ | # requirements at markets and supplies at factories. | ||
+ | # | ||
+ | # - Original formulation: | ||
+ | # Princeton University Press, Princeton, New Jersey, 1963. | ||
+ | # - Gams implementation: This formulation is described in detail in: | ||
+ | # Rosenthal, R E, Chapter 2: A GAMS Tutorial. In GAMS: A User's Guide. | ||
+ | # The Scientific Press, Redwood City, California, 1988. | ||
+ | # - JuMP implementation: | ||
- | """ | + | using CSV, DataFrames, GLPK, JuMP |
- | Basic example of transport model from GAMS model library translated to Pyomo | + | |
- | To run this you need pyomo and a linear solver installed. | + | # Sets |
- | When these dependencies are installed you can solve this example in one of | + | plants |
- | this ways (glpk is the default solver): | + | markets = [" |
- | ./ | + | # Parameters |
- | | + | a = Dict( # capacity of plant i in cases |
- | pyomo solve transport.py | + | " |
- | | + | " |
+ | ) | ||
+ | b = Dict( # demand at market j in cases | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | ) | ||
- | To display the results: | + | # distance in thousands of miles |
+ | d_table = CSV.read(IOBuffer(""" | ||
+ | plants | ||
+ | seattle | ||
+ | san_diego | ||
+ | """ | ||
+ | d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), | ||
- | cat results.json | + | f = 90 # freight in dollars per case per thousand miles |
- | cat results.yml (if PyYAML is installed on your system) | + | |
- | GAMS equivalent code is inserted as single-dash comments. The original GAMS code | + | c = Dict() # transport cost in thousands |
- | needs slighly different ordering | + | [ c[p,m] = f * d[p, |
- | http:// | + | |
- | Original problem formulation: | + | # Model declaration |
- | | + | trmodel = Model(GLPK.Optimizer) # transport model |
- | Princeton University Press, Princeton, New Jersey, 1963. | + | |
- | GAMS implementation: | + | |
- | Rosenthal, R E, Chapter 2: A GAMS Tutorial. In GAMS: A User's Guide. | + | |
- | The Scientific Press, Redwood City, California, 1988. | + | |
- | Pyomo translation: | + | |
- | Antonello Lobianco | + | |
- | This file is in the Public Domain | + | # Variables |
- | """ | + | @variables trmodel begin |
+ | x[p in plants, m in markets] >= 0 # shipment quantities | ||
+ | end | ||
- | # Import | + | # Constraints |
- | from pyomo.environ import * | + | @constraints trmodel begin |
+ | supply[p in plants], | ||
+ | sum(x[p,m] for m in markets) | ||
+ | demand[m in markets], | ||
+ | sum(x[p,m] for p in plants) | ||
+ | end | ||
- | # Creation of a Concrete Model | + | # Objective |
- | model = ConcreteModel() | + | @objective trmodel Min begin |
+ | sum(c[p, | ||
+ | end | ||
- | ## Define sets ## | + | print(trmodel) |
- | # Sets | + | |
- | # | + | |
- | # | + | |
- | model.i = Set(initialize=[' | + | |
- | model.j = Set(initialize=[' | + | |
- | ## Define parameters ## | + | optimize!(trmodel) |
- | # | + | status = termination_status(trmodel) |
- | # a(i) | + | |
- | # / | + | if status == MOI.OPTIMAL |
- | # san-diego | + | |
- | # b(j) | + | |
- | # / | + | |
- | # chicago | + | |
- | # topeka | + | |
- | model.a | + | |
- | model.b = Param(model.j, initialize={' | + | |
- | # Table d(i,j) | + | |
- | # new-york | + | |
- | # seattle | + | |
- | # san-diego | + | |
- | dtab = { | + | |
- | (' | + | |
- | (' | + | |
- | | + | |
- | (' | + | |
- | (' | + | |
- | (' | + | |
- | | + | |
- | model.d = Param(model.i, model.j, initialize=dtab, | + | |
- | # Scalar f freight in dollars per case per thousand miles /90/ ; | + | |
- | model.f = Param(initialize=90, doc=' | + | |
- | # Parameter c(i,j) transport cost in thousands of dollars per case ; | + | |
- | # c(i,j) = f * d(i,j) / 1000 ; | + | |
- | def c_init(model, | + | |
- | return model.f * model.d[i,j] / 1000 | + | |
- | model.c = Param(model.i, | + | |
- | ## Define variables ## | + | else |
- | # Variables | + | |
- | # x(i,j) shipment quantities in cases | + | |
- | # | + | end |
- | # Positive Variable x ; | + | |
- | model.x = Var(model.i, | + | |
- | + | ||
- | ## Define contrains ## | + | |
- | # supply(i) observe supply limit at plant i | + | |
- | # supply(i) .. sum (j, x(i,j)) =l= a(i) | + | |
- | def supply_rule(model, | + | |
- | return sum(model.x[i, | + | |
- | model.supply = Constraint(model.i, | + | |
- | # demand(j) | + | |
- | # demand(j) .. sum(i, x(i,j)) =g= b(j); | + | |
- | def demand_rule(model, | + | |
- | return sum(model.x[i, | + | |
- | model.demand = Constraint(model.j, | + | |
- | ## Define Objective and solve ## | ||
- | # cost define objective function | ||
- | # cost .. z =e= sum((i,j), c(i, | ||
- | # Model transport /all/ ; | ||
- | # Solve transport using lp minimizing z ; | ||
- | def objective_rule(model): | ||
- | return sum(model.c[i, | ||
- | model.objective = Objective(rule=objective_rule, | ||
- | |||
- | |||
- | ## Display of the output ## | ||
- | # Display x.l, x.m ; | ||
- | def pyomo_postprocess(options=None, | ||
- | model.x.display() | ||
- | |||
- | # This is an optional code path that allows the script to be run outside of | ||
- | # pyomo command-line. | ||
- | if __name__ == ' | ||
- | |||
- | #This replicates what the pyomo command-line tools does | ||
- | from pyomo.opt import SolverFactory | ||
- | import pyomo.environ | ||
- | opt = SolverFactory(" | ||
- | instance = model.create() | ||
- | results = opt.solve(instance) | ||
- | #sends results to stdout | ||
- | results.write() | ||
- | pyomo_postprocess(None, | ||
- | |||
# Expected result: | # Expected result: | ||
# obj= 153.675 | # obj= 153.675 | ||
Line 470: | Line 333: | ||
# | # | ||
</ | </ | ||
- | |||
- | |||
- | |||
- | ~~DISCUSSION~~ | ||
- | |||
- | |||
- | |||
~~DISCUSSION~~ | ~~DISCUSSION~~ | ||