Differences

This shows you the differences between two versions of the page.

Link to this comparison view

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 [2018/06/18 15:11] (current)
Line 35: Line 35:
 Pkg.add("​GLPKMathProgInterface"​) ​   # A lineaqr and MIP solver Pkg.add("​GLPKMathProgInterface"​) ​   # A lineaqr and MIP solver
 Pkg.add("​Ipopt"​) ​                   # A non-linear solver Pkg.add("​Ipopt"​) ​                   # A non-linear solver
-Pkg.add("​DataFrames"​) ​              # A library to deal with dataframes (R-like tabular data)+Pkg.add("​DataFrames"​) ​              # A library to deal with dataframes (R like tabular data)
 </​code>​ </​code>​
  
Line 44: Line 44:
 You will need to import as a minima the ''​JuMP''​ module. If you wish to specify a solver engine rather than letting JuMP select a suitable one, you will need to import also the module relative to the solver, e.g. ''​Ipopt''​ or  ''​GLPKMathProgInterface''​ You will need to import as a minima the ''​JuMP''​ module. If you wish to specify a solver engine rather than letting JuMP select a suitable one, you will need to import also the module relative to the solver, e.g. ''​Ipopt''​ or  ''​GLPKMathProgInterface''​
  
-<code julia> +<​code ​ julia> 
-# Import of the JuMP and DataFrames modules (the latter one just to import the data from a header-based table, as in the original trasnport example in GAMS +# Import of the JuMP and DataFrames modules (the latter one just to import the data from a header based table, as in the original trasnport example in GAMS 
 using JuMP, DataFrames using JuMP, DataFrames
 </​code>​ </​code>​
Line 56: Line 56:
    
 <code julia> <code julia>
-# Sets+Define sets # 
 +#  ​Sets 
 +#       ​i ​  ​canning plants ​  / seattle, san-diego / 
 +#       ​j ​  ​markets ​         / new-york, chicago, topeka / ;
 plants ​ = ["​seattle","​san_diego"​] ​         # canning plants plants ​ = ["​seattle","​san_diego"​] ​         # canning plants
 markets = ["​new_york","​chicago","​topeka"​] ​ # markets markets = ["​new_york","​chicago","​topeka"​] ​ # markets
Line 64: Line 67:
 ==== Definition of the "​parameters"​ ==== ==== Definition of the "​parameters"​ ====
  
-Capacity of plants and demand of markets are directly defined as dictionaries,​ while the distance is first read from a white-space separated table and then it is converted in a "​(plant,​ market) => value" dictionary.+Capacity of plants and demand of markets are directly defined as dictionaries,​ while the distance is first read as a DataFrame ​from a white-space separated table and then it is converted in a "​(plant,​ market) => value" dictionary.
  
 <code julia> <code julia>
-# Parameters+Define parameters # 
 +#   Parameters 
 +#       ​a(i) ​ capacity of plant i in cases 
 +#         / ​   seattle ​    350 
 +#              san-diego ​  ​600 ​ /
 a = Dict(              # capacity of plant i in cases a = Dict(              # capacity of plant i in cases
   "​seattle" ​  => 350,   "​seattle" ​  => 350,
   "​san_diego"​ => 600,   "​san_diego"​ => 600,
 ) )
 +
 +#       ​b(j) ​ demand at market j in cases
 +#         / ​   new-york ​   325
 +#              chicago ​    300
 +#              topeka ​     275  / ;
 b = Dict(              # demand at market j in cases b = Dict(              # demand at market j in cases
   "​new_york" ​ => 325,   "​new_york" ​ => 325,
Line 78: Line 90:
 ) )
  
-#  distance in thousands of miles+Table d(i,​j) ​ ​distance in thousands of miles 
 +#                    new-york ​      ​chicago ​     topeka 
 +#      seattle ​         2.5           ​1.7 ​         1.8 
 +#      san-diego ​       2.5           ​1.8 ​         1.4  ;
 d_table = wsv"""​ d_table = wsv"""​
 plants ​    ​new_york ​ chicago ​ topeka plants ​    ​new_york ​ chicago ​ topeka
Line 85: Line 100:
 """​ """​
 d = Dict( (r[:​plants],​m) => r[Symbol(m)] for r in eachrow(d_table),​ m in markets) d = Dict( (r[:​plants],​m) => r[Symbol(m)] for r in eachrow(d_table),​ m in markets)
 +# Here we are converting the table in a "​(plant,​ market) => distance"​ dictionary
 # r[:​plants]: ​  the first key, row field using a fixed header # r[:​plants]: ​  the first key, row field using a fixed header
 # m:            the second key # m:            the second key
 # r[Symbol(m)]:​ the value, the row field with a dynamic header # r[Symbol(m)]:​ the value, the row field with a dynamic header
  
 +# 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,​j) ​ transport cost in thousands of dollars per case ;
 +#            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] ​
 </​code>​ </​code>​
 +The above code take advantage of [[http://​docs.julialang.org/​en/​stable/​manual/​arrays/#​comprehensions|List Comprehensions]],​ a powerful feature of the Julia language that provides a concise way to loop over a list.
 +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[:​plants],​m) => r[Symbol(m)]
 +  end
 +end
 +</​code>​
 +Using List Comprehension is however quicker to code and more readable.
  
  
Line 105: Line 135:
  
 <code julia> <code julia>
-# Model declaration +# Model declaration ​(transport model) 
-trmodel = Model() ​# transport model+trmodel = Model() ​ 
 </​code>​ </​code>​
  
 ==== 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 ## 
 +#  ​Variables 
 +#       ​x(i,​j) ​ shipment quantities in cases 
 +#       ​z ​      total transportation costs in thousands of dollars ; 
 +#  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 161:
  
 <code julia> <code julia>
-Constraints+## Define contrains ## 
 +# supply(i) ​  ​observe supply limit at plant i 
 +# supply(i) .. sum (j, x(i,j)) =l= a(i) 
 +# demand(j) ​  ​satisfy demand at market j ;   
 +# demand(j) .. sum(i, x(i,j)) =g= b(j);
 @constraints trmodel begin @constraints trmodel begin
     supply[p in plants], ​  # observe supply limit at plant p     supply[p in plants], ​  # observe supply limit at plant p
Line 151: Line 191:
 <code julia> <code julia>
 print(trmodel) print(trmodel)
-</code+</code>
  
 ==== Resolution of the model ==== ==== Resolution of the model ====
Line 159: Line 199:
 <code julia> <code julia>
 status = solve(trmodel) status = solve(trmodel)
-</julia>+</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 ''​getvalue(var_name)'',​ you can just ''​println(getvalue(x))''​ to get a basic output.\\
 +Notice that you can also easily retrieve the dual value associated to the constraint with ''​getdual(constraint_name)''​.
  
 <code julia> <code julia>
Line 167: Line 209:
     println("​Objective value: ", getobjectivevalue(trmodel))     println("​Objective value: ", getobjectivevalue(trmodel))
     println(getvalue(x))     println(getvalue(x))
 +    println("​Shadow prices of supply:"​)
 +    [println("​$p = $(getdual(supply[p]))"​) for p in plants]
 +    println("​Shadow prices of demand:"​)
 +    [println("​$m = $(getdual(demand[m]))"​) for m in markets]
 else else
     println("​Model didn't solved"​)     println("​Model didn't solved"​)
Line 173: Line 219:
 </​code>​ </​code>​
  
-==== 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,​ we are passing to the constructor the initial elements of the set and a documentation string to keep track on what our set represents: 
-<code python> 
-## Define sets ## 
-#  Sets 
-#       ​i ​  ​canning plants ​  / seattle, san-diego / 
-#       ​j ​  ​markets ​         / new-york, chicago, topeka / ; 
-model.i = Set(initialize=['​seattle','​san-diego'​],​ doc='​Canning plans'​) 
-model.j = Set(initialize=['​new-york','​chicago',​ '​topeka'​],​ doc='​Markets'​) 
-</​code>​ 
  
-==== Parameters ​==== +==== Editing and running the script ​==== 
-Parameter objects are created specifying the sets over which they are defined and are initialised with either a python dictionary or a scalar: +Differently from GAMS you can use whatever editor environment you wish to code a JuMP scriptIf you don't need debugging features, ​simple text editor like Notepad++ ​(in windows), gedit or kate (in Linuxwill sufficeThey already have syntax highlight for Julia.\\ 
-<code python>​ +If you want advanced features and debugging capabilities you can use a dedicated Julia IDElike e.g[[http://junolab.org/|Juno]].
-## Define parameters ## +
-#   ​Parameters +
-#       a(i)  capacity of plant i in cases +
-#         / ​   seattle ​    350 +
-#              san-diego ​  ​600 ​ / +
-#       ​b(j) ​ demand at market j in cases +
-#         / ​   new-york ​   325 +
-#              chicago ​    300 +
-#              topeka ​     275  / ; +
-model.a = Param(model.i, initialize={'​seattle':​350,'​san-diego':​600},​ doc='​Capacity of plant i in cases') +
-model.b = Param(model.jinitialize={'​new-york':​325,'​chicago':​300,'​topeka':​275},​ doc='​Demand at market j in cases'​) +
-#  Table d(i,j)  distance ​in thousands of miles +
-#                    new-york ​      ​chicago ​     topeka +
-#      seattle ​         2.5           ​1.7 ​         1.8 +
-#      san-diego ​       2.5           ​1.8 ​         1.4  ; +
-dtab = { +
-    ('​seattle', ​ '​new-york'​: 2.5, +
-    ('​seattle', ​ '​chicago'​) ​ : 1.7, +
-    ​('​seattle'​ '​topeka'​) ​  : 1.8, +
-    ('​san-diego','​new-york'​) : 2.5, +
-    ('​san-diego','​chicago'​)  ​1.8, +
-    ('​san-diego','​topeka'​) ​  : 1.4, +
-    } +
-model.d = Param(model.i,​ model.j, initialize=dtab,​ doc='​Distance in thousands of miles'​) +
-#  Scalar f  freight in dollars per case per thousand miles  ​/90+
-model.f = Param(initialize=90,​ doc='​Freight in dollars per case per thousand miles'​) +
-</code> +
-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,​j) ​ transport cost in thousands of dollars per case ; +
-#            c(i,j) = f * d(i,j) / 1000 ; +
-def c_init(model,​ i, j): +
-  return model.f * model.d[i,j/ 1000 +
-model.c = Param(model.i,​ model.j, initialize=c_init,​ doc='​Transport cost in thousands of dollar per case'​) +
-</​code>​+
  
-==== Variables ==== +If you are using instead the Julia console ​you ​can run the script as ''julia transport.jl''​.
-Similar to parametersvariables 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 +
-#       ​x(i,​j) ​ shipment quantities in cases +
-#       ​z ​      total transportation costs in thousands of dollars ; +
-#  Positive Variable x ; +
-model.x = Var(model.i,​ model.j, bounds=(0.0,​None),​ doc='Shipment quantities in case'+
-</​code>​+
  
-==== Constrains ​==== +===== Further help ===== 
-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. ​  +Documentation ​of JuMP is available from [[https://jump.readthedocs.io/en/latest/|this page]]. However if you want to do serious things with juMPit is most likely that you will have to either look at the source ​code or consult the [[https://discourse.julialang.org/​c/​domain/​opt|discussion forum]].
-<code python>​ +
-## Define contrains ## +
-# supply(i) ​  ​observe supply limit at plant i +
-# supply(i) .. sum (j, x(i,​j)) ​=la(i) +
-def supply_rule(model,​ i): +
-  return sum(model.x[i,​j] for j in model.j) <= model.a[i] +
-model.supply = Constraint(model.i,​ rule=supply_rule,​ doc='​Observe supply limit at plant i') +
-# demand(j) ​  ​satisfy demand at market j ;   +
-# demand(j) .. sum(i, x(i,j)) =g= b(j); +
-def demand_rule(model,​ j): +
-  return sum(model.x[i,​j] for i in model.i) >= model.b[j] ​  +
-model.demand = Constraint(model.j,​ rule=demand_rule,​ doc='​Satisfy demand at market j') +
-</​code>​ +
-The above code take advantage ​of [[https://docs.python.org/2/tutorial/datastructures.html#​list-comprehensions|List Comprehensions]], a powerful feature of the python language that provides a concise way to loop over a list. +
-If we take the supply_rule as examplethis 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,​ i): +
-  supply = 0.+
-  for j in model.j: +
-    supply += model.x[i,j] +
-  return supply <= model.a[i] +
-</​code>​ +
-Using List Comprehension is however quicker to code and more readable.+
  
-==== Objective & solving ==== +Happy modelling with JuMP ;-)
-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,​j)*x(i,​j)) ​; +
-#  Model transport /all/ ; +
-#  Solve transport using lp minimizing z ; +
-def objective_rule(model):​ +
-  return sum(model.c[i,​j]*model.x[i,​j] for i in model.i for j in model.j) +
-model.objective = Objective(rule=objective_rule,​ sense=minimize,​ doc='​Define objective function'​) +
-</​code>​ +
-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,​kj]*model.x[ki,​kj] +
-  return obj +
-</​code>​+
  
-==== Retrieving the output ​==== +===== Complete script ​=====
-To retrieve the output and do something with it (either to just display it -like we do here-, to plot a graph with [[http://​matplotlib.org|matplotlib]] or to save it in a csv file) we use the ''​pyomo_postprocess()''​ function.\\ +
-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, instance=None, results=None):​ +
-  model.x.display() +
-</​code>​ +
-We can print model structure information with ''​model.pprint()''​ ("​pprint"​ stand for "​pretty print"​).\\ +
-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 ​==== +Here is the complete ​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.\\ +
-If you want advanced features and debugging capabilities you can use a dedicated Python IDE, like e.g. [[https://​code.google.com/​p/​spyderlib/​|Spyder]].+
  
-You will normally run the script as ''​pyomo solve --solver=glpk transport.py''​. You can output solver specific output adding the option ''​--stream-output''​.\\ +<code Julia> 
-If you want to run the script as ''​python ​transport.py''​ add the following lines at the end:\\ +# Transposition in JuMP of the basic transport ​model used in the GAMS tutorial 
-<code python> + 
-# This is an optional code path that allows the script to be run outside of +# This problem finds a least cost shipping schedule ​that meets 
-pyomo command-line For example: ​ python transport.py +requirements at markets and supplies at factories
-if __name__ == '​__main__':​ + 
-    #This replicates what the pyomo command-line tools does +# - Original formulation:​ Dantzig, G B, Chapter 3.3In Linear Programming and Extensions
-    from pyomo.opt import SolverFactory +# Princeton University Press, Princeton, New Jersey, 1963
-    import pyomo.environ +- Gams implementation:​ This formulation is described in detail in: 
-    opt = SolverFactory("​glpk"​) +# Rosenthal, R E, Chapter 2: A GAMS Tutorial. In GAMS: A User's Guide
-    instance = model.create() +# The Scientific PressRedwood CityCalifornia, 1988. 
-    ​results = opt.solve(instance) +# - JuMP implementation:​ Antonello Lobianco
-    #sends results to stdout +
-    ​results.write() +
-    ​pyomo_postprocess(Noneinstanceresults) +
-</​code>​+
  
-Finallyif you are very lazy and want to run the script with just ./​transport.py (and you are in Linux) add the following lines at the top: +using JuMPDataFrames
-<code python>​ +
-#​!/​usr/​bin/​env python +
-# -*- coding: utf-8 -*- +
-</​code>​+
  
-===== Further help ===== +# Sets 
-Documentation of pyomo is available from [[https://​software.sandia.gov/​trac/​coopr/​wiki/​Documentation|this page]]. However if you want to do serious things with pyomo, it is most likely that you will have to either look at the source code or consult the [[https://​groups.google.com/​forum/​#!forum/​coopr-forum|mailing list]].+plants  ​["​seattle","​san_diego"​] ​         # canning plants 
 +markets ​= ["​new_york","​chicago","​topeka"​ markets
  
-Happy modelling with pyomo ;-) +# Parameters 
-===== Complete script =====+Dict(              # capacity of plant i in cases 
 +  "​seattle" ​  => 350, 
 +  "​san_diego" ​=> 600, 
 +
 +Dict(              # demand at market j in cases 
 +  "​new_york"  ​=> 325, 
 +  "​chicago" ​  => 300, 
 +  "​topeka" ​   ​=> 275, 
 +)
  
-Here is the complete script: ​ +#  distance in thousands of miles 
- +d_table = wsv"""​ 
-<code python> +plants ​    ​new_york ​ chicago ​ topeka 
-#​!/​usr/​bin/​env python +seattle ​   2.5       ​1.7 ​     1.8 
-# -*- coding: utf--*- +san_diego ​ 2.5       1.     1.4
- +
 """​ """​
-Basic example of transport model from GAMS model library translated to Pyomo +d = Dict( (r[:plants],m) => r[Symbol(m)] for in eachrow(d_table), in markets)
-  +
-To run this you need pyomo and a linear solver installed. +
-When these dependencies are installed you can solve this example in one of +
-this ways (glpk is the default solver): +
-  +
-    ./​transport.py ​(Linux only) +
-    python transport.py +
-    pyomo solve transport.py +
-    pyomo solve --solver=glpk transport.py +
-  +
-To display the results: +
-  +
-    cat results.json +
-    cat results.yml (if PyYAML is installed on your system) +
-  +
-GAMS equivalent code is inserted as single-dash comments. The original GAMS code +
-needs slighly different ordering of the commands and it's available at +
-http://​www.gams.com/​mccarl/​trnsport.gms +
-  +
-Original problem formulation:​ +
-    Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions. +
-    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 +
-"""​ +
-  +
-# Import +
-from pyomo.environ import * +
-  +
-# Creation of a Concrete Model +
-model = ConcreteModel() +
-  +
-## Define sets ## +
-#  Sets +
-#       ​i ​  ​canning ​plants ​  / seattle, san-diego / +
-#       ​j ​  ​markets ​         / new-york, chicago, topeka / ; +
-model.i = Set(initialize=['​seattle','​san-diego'​], doc='​Canning plans') +
-model.j = Set(initialize=['​new-york','​chicago',​ '​topeka'​],​ doc='​Markets'​) +
-  +
-## Define parameters ## +
-#   ​Parameters +
-#       a(i ​capacity of plant i in cases +
-#         / ​   seattle ​    350 +
-#              san-diego ​  ​600 ​ / +
-#       ​b(j) ​ demand at market j in cases +
-#         / ​   new-york ​   325 +
-#              chicago ​    300 +
-#              topeka ​     275  / ; +
-model.a = Param(model.i,​ initialize={'​seattle':​350,'​san-diego':​600},​ doc='​Capacity of plant i in cases'​) +
-model.b = Param(model.j,​ initialize={'​new-york':​325,'​chicago':​300,'​topeka':​275},​ doc='​Demand at market j in cases'​) +
-#  Table d(i,​j) ​ distance in thousands of miles +
-#                    new-york ​      ​chicago ​     topeka +
-#      seattle ​         2.5           ​1.7 ​         1.8 +
-#      san-diego ​       2.5           ​1.8 ​         1.4  ; +
-dtab = { +
-    ('​seattle', ​ '​new-york'​) : 2.5, +
-    ('​seattle', ​ '​chicago'​) ​ : 1.7, +
-    ('​seattle', ​ '​topeka'​) ​  : 1.8, +
-    ('​san-diego','​new-york'​) : 2.5, +
-    ('​san-diego','​chicago'​) ​ : 1.8, +
-    ('​san-diego','​topeka'​) ​  : 1.4, +
-    } +
-model.d = Param(model.i,​ model.j, initialize=dtab,​ doc='​Distance in thousands of miles'​) +
-#  Scalar f  freight in dollars per case per thousand miles  /90/ ; +
-model.f = Param(initialize=90,​ doc='​Freight in dollars per case per thousand miles'​) +
-#  Parameter c(i,​j) ​ transport cost in thousands of dollars per case ; +
-#            c(i,j) = f * d(i,j) / 1000 ; +
-def c_init(model,​ i, j): +
-  return model.f * model.d[i,​j] / 1000 +
-model.c = Param(model.i,​ model.j, initialize=c_init,​ doc='​Transport cost in thousands of dollar per case'​) +
-  +
-## Define variables ## +
-#  Variables +
-#       ​x(i,​j) ​ shipment quantities in cases +
-#       ​z ​      total transportation costs in thousands of dollars ; +
-#  Positive Variable x ; +
-model.x = Var(model.i,​ model.j, bounds=(0.0,​None),​ doc='​Shipment quantities in case'​) +
-  +
-## Define contrains ## +
-# supply(i) ​  ​observe supply limit at plant i +
-# supply(i) .. sum (j, x(i,j)) =l= a(i) +
-def supply_rule(model,​ i): +
-  return sum(model.x[i,​j] for in model.j) <= model.a[i] +
-model.supply = Constraint(model.i, rule=supply_rule,​ doc='​Observe supply limit at plant i') +
-# demand(j) ​  ​satisfy demand at market j ;   +
-# demand(j) .. sum(ix(i,j)) =g= b(j); +
-def demand_rule(model,​ j): +
-  return sum(model.x[i,​j] for i in model.i) >= model.b[j] ​  +
-model.demand = Constraint(model.j,​ rule=demand_rule,​ doc='​Satisfy demand at market j')+
  
-## Define Objective and solve ## +f = 90 freight in dollars per case per thousand miles 
-#  cost        define objective function + 
-#  cost ..        z  =e sum((i,j), c(i,​j)*x(i,​j)) ; +Dict() # transport ​cost in thousands of dollars per case 
- ​Model ​transport ​/all/ +c[p,m= f d[p,m/ 1000 for in plants, m in markets]
-#  Solve transport using lp minimizing z ; +
-def objective_rule(model):​ +
-  return sum(model.c[i,j]*model.x[i,j] for in model.i for j in model.j) +
-model.objective = Objective(rule=objective_rule,​ sense=minimize,​ doc='​Define objective function'​)+
  
-  +# Model declaration 
-## Display of the output ​## +trmodel = Model() ​transport model 
-Display ​x.lx.; + 
-def pyomo_postprocess(options=Noneinstance=Noneresults=None): +Variables 
-  ​model.x.display() +@variables trmodel begin 
-  +    x[p in plants, m in markets] >= 0 shipment quantities in cases 
-# This is an optional code path that allows the script to be run outside of +end 
-# pyomo command-line. ​ For example: ​ python transport.py + 
-if __name__ ​== '​__main__'​: +Constraints 
-  +@constraints trmodel begin 
-    ​#This replicates what the pyomo command-line tools does +    supply[p in plants], ​  observe supply limit at plant p 
-    ​from pyomo.opt import SolverFactory +        sum(x[p,m] for m in markets) ​ <=  a[p] 
-    ​import pyomo.environ +    ​demand[m in markets], ​ # satisfy demand at market m 
-    ​opt = SolverFactory("glpk") +        sum(x[p,m] for p in plants) ​ > ​b[m] 
-    ​instance ​model.create() +end 
-    ​results = opt.solve(instance+ 
-    ​#sends results to stdout +# Objective 
-    ​results.write() +@objective trmodel Min begin 
-    ​pyomo_postprocess(None, instance, results)+    sum(c[p,m]*x[p,m] for p in plants, m in markets
 +end 
 + 
 +print(trmodel
 + 
 +status = solve(trmodel) 
 + 
 +if status ​== :Optimal 
 +    ​println("​Objective value: ", getobjectivevalue(trmodel)) 
 +    ​println("​Shipped quantities: ") 
 +    ​println(getvalue(x)) 
 +    ​println("Shadow prices of supply:") 
 +    ​[println("​$p ​$(getdual(supply[p]))"​for p in plants] 
 +    ​println("​Shadow prices of demand:"​
 +    ​[println("​$m = $(getdual(demand[m]))"​) for m in markets] 
 + 
 +else 
 +    ​println("Model didn't solved"​
 +    ​println(status) 
 +end
    
 # Expected result: # Expected result:
Line 470: Line 326:
 #​['​san-diego','​topeka'​] ​  = 275 #​['​san-diego','​topeka'​] ​  = 275
 </​code>​ </​code>​
- 
- 
- 
-~~DISCUSSION~~ 
- 
- 
- 
  
 ~~DISCUSSION~~ ~~DISCUSSION~~
  
personal/blog/2017/0203_jump_for_gams_users.1486129506.txt.gz · Last modified: 2018/06/18 15:10 (external edit)
CC Attribution-Noncommercial-Share Alike 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0