Modeling OR Problems in SAGE

Sage is an open source project providing a comprehensive environment for mathematical computing. It's a big package with big ambitions "Creating a viable free open source alternative to Magma, Maple, Mathematica and Matlab."  See SIAM Review (Vol 41, No. 4, 2009) for a  recent review.

Sage incorporates Python as a scripting language from which one can generate objects MixedIntegerLinearProgram objects.  Both GLPK and CBC from the COIN-OR project can be installed as optional Sage packages and used as solvers.  Coupled with the native features of Python, this approaches provides an interesting alternative to developing models using AMPL/GMPL.

# Dictionary of task durations indexed by (Job, Machine) tuples
dur = {
    ('Paper 1','Blue')   : 45,
    ('Paper 1','Yellow') : 10,
    ('Paper 2','Blue')   : 20,
    ('Paper 2','Green')  : 10,
    ('Paper 2','Yellow') : 34,
    ('Paper 3','Blue')   : 12,
    ('Paper 3','Green')  : 17,
    ('Paper 3','Yellow') : 28 }

# Task sequencing indicated by pairs of (Job, Machine) tuples
order = [
    (('Paper 1','Blue'),  ('Paper 1','Yellow')),
    (('Paper 2','Green'), ('Paper 2','Blue')),
    (('Paper 2','Blue'),  ('Paper 2','Yellow')),
    (('Paper 3','Yellow'),('Paper 3','Blue')),
    (('Paper 3','Blue'),  ('Paper 3','Green')) ]

# TASKS are a list of (job,machine) tuples found from the keys for dur{}
TASKS = dur.keys()

# JOBS and MACHS are unduplicated lists of jobs and machines 
JOBS = sorted(list(set(zip(*TASKS)[0])))
MACHS = sorted(list(set(zip(*TASKS)[1])))

# Create MILP
JobShop = MixedIntegerLinearProgram(maximization=False)

# The objective is to minimize Makespan
MakeSpan = JobShop.new_variable()
JobShop.set_objective(MakeSpan[0])

# MakeSpan is an upper bound on all tasks
x = JobShop.new_variable()
[JobShop.add_constraint(x[t] + dur[t] - MakeSpan[0], max = 0) for t in TASKS]

# Tasks must be complete in the specified order 
[JobShop.add_constraint(x[s] + dur[s] - x[t], max=0) for s,t in order]

# Disjunctive constraints to avoid machine conflicts. This uses a
# 'Big M' technique and a binary variable where y[s][t] = 0 implies
# task s must finish before task t can start and y[s][t] = 1 implies
# task t finishes before task s starts.

CONFLICTS = [(s,t) for s in TASKS for t in TASKS if s<t and s[1]==t[1]]

y = JobShop.new_variable(dim=2,vtype=JobShop.__BINARY)

BigM = sum(dur[t] for t in TASKS)

[JobShop.add_constraint(BigM*y[s][t] + x[t] - x[s] - dur[s], min=0)
    for s,t in CONFLICTS]
[JobShop.add_constraint(BigM*(1-y[s][t]) + x[s] - x[t] - dur[t], min=0)
    for s,t in CONFLICTS]
JobShop.solve('Coin')

Gambling with Stochastic Dynamic Programming

Stochastic dynamic programming is a technique for multi-period decision making under uncertainty. A classic illustration is the case of a risk neutral gambler entering a game with a stake x and a goal of leaving the game with a stake N > x.  The gambler places a wager from his or her current stake, then either wins an equal amount with probability p or loses the wager with probability q. The game continues until the gambler either reaches the goal or goes bust. The gambler cares only about the expected value and not about the range of possible outcomes and therefore is said to be risk neutral.  The problem is to find a betting strategy maximizing the expected value of the game.

The usual setup for solving this problem involves Bellman's 'principle of optimality'. The optimality criterion establishes a recursive expression for the expected value of the game that is a function of the stake. The 'Bellman equation' can be solved for an optimal strategy by a well-known technique called 'policy iteration'.  Alternatively, the same equation can be solved by linear programming as demonstrated in the GMPL file RNGambling.mod.

A related example is the case of the risk averse gambler.  The risk averse gambler enters the game for a predetermined number of rounds with a goal of maximizing the expected 'utility' of the final outcome.  Utility is a concave function of wealth that more heavily penalizes losses than gains. The GMPL model RAGambling.mod demonstrates that this problem can also be solved with linear programming. This example is a closely related to the well known 'Kelly Criterion' (or here for a popular account).

Example of Business Planning through Stochastic Programming

The final project/exam for the course concerned a business planning problem that could be formulated as a two stage stochastic program.  I've posted here a variation of that problem as an additional example of GMPL modeling: BusinessPlanning.mod.

Final Exam, Assignments, and Assignment 7

Given the collision of so many end-of-semester activities, let's forego any new assignments for this week. Assignment 7 will be due Wednesday followed by a take-home final that will be due noon on Thursday, December 17th.  We'll take more about the final during Wednesday's class session.

Regarding Assignment 7, please feel free to use the GMPL model for portfolio optimization using the MAD criterion.  You'll need to adapt it to the problem data given in the assignment, and in particular deal with the portfolio constraint in part (b) of the problem.

Portfolio Optimization with a Mean Absolute Deviation Criterion

A question came up after class on whether GLPK could be used for portfolio optimization.  At first glance the answer is no because the classic Markowitz formulation for portfolio optimization with a  mean-variance criterion leads to a quadratic program.  GLPK is a linear programming solver that doesn't generally handle quadratic programs. There are workarounds of limited utility, which is demonstrated in a post below.  But in general, GLPK is not suited to quadratic programming.

However, changing the optimization criterion leads to some interesting and useful tools for portfolio optimization that can be expressed as a linear programs. An example is a criterion based on the mean of the absolute deviation (MAD) of return as a risk measure, an idea attributed to Konno and Yamazaki (1991).  This approach also has some other interesting features, including the direct use of return data either from simulation or historical records. Solutions for the MAD criterion exhibit second order stochastic dominance - a a theoretically important feature for risk measures that is not shared by the Markowitz formulation.

To demonstrate the concept, I've prepared a sample GMPL file PortfolioMAD.mod that computes an optimal portfolio using the MAD criterion.  The calculation shows how to compute random samples from multivariable Normal distribution in GMPL.  This requires an implementation of a Cholesky factorization in GMPL which is given in the code.

Notes on Assignment 6: Australian Motors Case Study

I've been responding to questions regarding the Australian Motors (AM) case study and wanted to pass on a couple of observations.

First, the purpose of this exercise is to implement the business model outlined in Exhibit 8.  AM leases trucks on 2, 3, 4, and 5 month contracts. It makes money by renting those trucks to customers for short term needs.  AM is like a rental car company that obtains it cars through leasing.  Your task is to implement a business model that AM can use to schedule leases for maximum profit.

Second, the model described in Exhibit 8 assigns no cost to AM for leasing surplus vehicles for its fleet that don't, in turn, get rented out to customers. This unrealistic assumption leads to a trivial solution.  You'll get a more realistic solution if you modify the objective function to include an additional cost of $500 per month for each surplus vehicle in the inventory.

Take some time to look over your solution.  In particular, the marginal values on the inequality constraints provides much insight into this business case.  For example, which months and demand segments would your prioritize for a marketing effort intended to increased profits?  Do you see months and demand segments where increasing demand by one unit would increase profits by more than the monthly contribution of that unit?  Why does that occur?

Notes on Assignment 7: Portfolio Optimization in GMPL/GLPK

I was reminded that not all of you are experienced with Matlab, and was asked whether it would be possible to do the portfolio optimization problem in Assignment 7 using GMPL/GLPK.  The answer is a qualified yes.

The issue is that the Mean-Variance (MV) formulation for portfolio optimization leads to a quadratic programming problem.  GLPK is for linear programs, thus the quadratic objective cannot be directly introduced to GLPK.  For the first part of this assignment, however, you can work around the issue by using the first order necessary conditions that define a solution for the MV portfolio.

Look at page 30 of the notes for Lecture 7 entitled "General Case". There you will find linear equations which constrain solutions to the MV problem.  Two of these are simple equality constraints that can be directly translated into GMPL using portfolio weights as the decision variables.  The remaining n linear equations involves the weights and two additional decision variables for the Lagrange multipliers.  Adding those two decision variables to the model allows you to translate these equations into a system of GMPL constraints.  You won't need an objective function because the linear constraints are sufficient to determine a unique solution for the MV portfolio.

For your convenience, I've put a GMPL data file called PortfolioMVO.dat on the downloads list with the necessary data for Assignment 7.

The trouble you'll face is that the second part of Assignment introduces an additional inequality constraint. This change requires you to go back and establish KKT necessary conditions that include this new constraint.  In the end you will have have a somewhat more complex GMPL model with additional decision variables and some discrete logic.  If you're interested in doing this, let me know and I'd be glad to work you through the mechanics.