Coopr: A Python Repository for Optimization

1. Coopr: A Python Repository for Optimization William E. Hart Sandia National Laboratories wehart@sandia.gov

2. Overview Why Python? Python Optimization Resources Coopr Overview Pyomo Modeling Tool

3. Why Python? • Features • A rich set of datatypes, support for object oriented programming, namespaces, exceptions, dynamic loading, etc. • Includes a large number of useful modules. • Syntax • Python has a clean syntax: no weird symbols (e.g. $, %, @). • Natural support for procedural and object-oriented software design. • Extendability and Customization • Simple model for loading Python code developed by a user. • Can easily integrate libraries that optimize compute kernels • Python can dynamically integrate libraries

4. Why Python? • Interpreter • A powerful interactive interpreter that allows real-time code development and encourages experimentation with Python software. • Documentation • Extensive online documentation and many excellent Python books • Support and Stability • Python is highly stable, and well supported. • Portability • Python is available on a wide range of compute platforms. • Open Source License • Python is freely available, with a liberal open source license

5. Use of Python in Scientific Computing • NumPy • SciPy • SAGE • JPL • PMV • AutoDock • IPython

6. Python Optimization Overview • CVXOPT • SciPy • OpenOpt • NLPy • Pyipopt • PuLP • POAMS

7. Coopr Overview Coopr is composed of a variety of Python Packages: • Coopr.util • Utility library • Coopr.opt • Generic optimization interfaces • Coopr.pyomo • A mathematical programming modeling tool • Coopr.exact • A tool for managing computational experiments

8. Pyomo Idea: support mathematical modeling of integer programs in Python

9. Application Example Min sum {j in P} c_j X_j s.t. sum {j in P} X_j/a_j 0 <= X_j <= u_j

10. Mathematical Modeling Languages Examples: • AMPL/MathProg • GAMS • AIMMS • OptimJ • … Impact: • Use a simple, natural syntax to describe mathematical models • Can formulate large models with a concise syntax • Separate modeling and data declarations • Support sophisticated data indexing to make modeling structuring very flexible • Support data import and export in commonly used formats • Include tools for debugging a model

11. AMPL Model set P; param a {j in P}; param b ; param c {j in P}; param u {j in P}; var X {j in P}; maximize Total_Profit : sum {j in P} c[j] * X[j]; subject to Time : sum {j in P} (1/ a[j]) * X[j] <= b; subject to Limit {j in P}: 0 <= X[j] <= u[j];

12. Example of AMPL Data data; set P := bands coils; param: a c u := bands 200 25 6000 coils 140 30 4000 ; param b := 40;

13. Pyomo Motivation • Support for flexible programming language features • Open specification of semantics/syntax • Supports reuse (e.g. class definitions, importing common routines) • Support for multiple model instances • Support for expressive feedback from optimizers • Rely on mature Python infrastructure for common data processing • data feed and cleaning • graphical reporting • client-server interaction • web-enabling of applications

14. Pyomo Motivation • Python extensions enable performance optimization of the modeling process • Syntax of modeling languages tacitly limit the degree to which the modeling can be customized • Can easily interface with solvers • Integrate python solvers • Linkage to high-performance numerical libraries • Support mathematical modeling with an open-source framework • Encourages updates by users • Enables community-wide critique of modeling capabilities

15. Pyomo Model (1) # # Coopr import # from coopr.pyomo import ¤ # # Setup the model # model = Model() # # Declare sets, parameters and variables # model.P = Set() model.a = Param(model.P) model.b = Param() model.c = Param(model.P) model.u = Param(model.P) model.X = Var(model.P)

16. Pyomo Model (2) # # Declare objective rule and create objective object # def Objective_rule(model): ans = 0 for j in model.P: ans = ans + model.c[j] * model.X[j] return ans model.Total_Profit = Objective(rule=Objective_rule, sense=maximize)

17. Pyomo Model (3) # # Declare constraint rules and create objective objects # # Time # def Time_rule(model): ans = 0 for j in model.P: ans = ans + (model.X[j]/model.a[j]) return ans < model.b model.Time = Constraint(rule=Time_rule) # # Limit # def Limit_rule(j, model): return(0, model.X[j], model.u[j]) model.Limit = Constraint(model.P, rule=Limit_rule)

18. Solving AMPL Model % ampl ampl: model prod.mod; ampl: data prod.dat; option solver PICO; solve; No integer variables... solving the LP normally. LP value= 192000 CPU RunTime= 0 CPU TotalTime= 0.109 WallClock TotalTime= 0.109375 PICO Solver: final f = 192000.000000

19. Solving Pyomo Model (1) % prod.py prod.dat ========================================================== --- Solver Results --- ========================================================== ---------------------------------------------------------- ------ Problem Information ------ ---------------------------------------------------------- name: None num_constraints: 5 num_nonzeros: 6 num_objectives: 1 num_variables: 2 sense: maximize upper_bound: 192000

20. Solving Pyomo Model (2) ---------------------------------------------------------- ------ Solver Information ------ ---------------------------------------------------------- error_rc: 0 nbounded: None ncreated: None status: ok systime: None usrtime: None ---------------------------------------------------------- ------ Solution 0 ---------------------------------------------------------- gap: 0.0 status: optimal value: 192000 Primal Variables X_bands_ 6000 X_coils_ 1400 Dual Variables c_u_Limit_1 4 c_u_Time_0 4200 ----------------------------------------------------------

21. Solving a Pyomo Model within Python # # Imports import prod from coopr.pyomo import * # # Configure Coopr coopr.opt.config().configure() # # Create the model instance instance = prod.model.create("prod.dat") # # Setup the optimizer opt = solvers.SolverFactory("glpk") # # Optimize results = opt.solve(instance) # # Write the output results.write(num=1)

22. AMPL/Pyomo Comparison • Pyomo object/constraint declarations are more verbose • Typically requires the use of a temporary function • Pyomo declarations explicitly refer to models • Can declare multiple models and • Pyomo can apply solvers that do not recognize the *.nl format • Coopr.opt supports automatic conversion to solver-specific format • LP format, MPS format • Pyomo (currently) only supports linear models • Linear programs and mixed-integer linear programs • Pyomo does not include preprocessing of LP/MILP models

23. But is it easy to use??? Perhaps… Case Study: Cargo Logistics Model • Plan logistics for a fleet of cargo aircraft • Collaboration between Lockheed-Martin and Sandia Labs • Initial mathematical prototype developed between LM and Sandia • Model was refined and extended at LM • Internal discussions with another business unit • LM build the Pyomo model and debugged it • Yes, there were a few Pyomo bugs, but not too many… Observations: • A non-expert user was able to use Pyomo without difficulty • Excellent Python documentation used by user • Excellent Python unit testing tools key to the stability of Pyomo

24. Future Work • Other solver interfaces for coopr.opt • Direct solver interfaces • Interface with SciPy, OpenOpt, etc… • Integration with SAGE • E.g. use formulas defined in SAGE • Generation of nonlinear models in Pyomo • Support for additional modeling constructs in Pyomo • E.g. stochastic programs • Piecewise linear problems • Optimization visualization and analysis in coopr.opt • Integration with PuLP and POAMS