Solving Linear and Integer Programs Robert E. Bixby ILOG, Inc. - - PowerPoint PPT Presentation

Solving Linear and Integer Programs

Robert E. Bixby

ILOG, Inc. and Rice University

2

Outline

Linear Programming:

Introduction to basic LP, including duality Primal and dual simplex algorithms Computational progress in linear programming Implementing the dual simplex algorithm

Mixed-Integer Programming:

3

Some Basic Theory

4

Where c∈Rn, b∈Rm, A∈Rm×n, and x is a vector of n variables. cTx is known as the objective function, Ax=b as the constraints, and x ≥ 0 as the nonnegativity conditions. b is called the right-hand side. (P) Minimize cTx Subject to Ax = b x ≥ 0 A linear program (LP) in standard form is an

• ptimization problem of the form

Linear Program – Definition

5

In this context, (P) is referred to as the primal linear program.

Minimize cTx Subject to Ax = b x ≥ 0

(D) Maximize bTπ Subject to ATπ ≤ c π free The dual (or adjoint) linear program corresponding to (P) is the optimization problem

Primal

Dual Linear Program – Definition

6

If bTπ = cTx, then x is optimal for (P) and π is

• ptimal for (D); moreover, if either (P) or (D) is

unbounded, then the other problem is infeasible. bTπ ≤ cTx Let x be feasible for (P) and π feasible for (D). Then

Weak Duality Theorem

(von Neumann 1947)

Proof: πTAx πTb =

Ax = b

≤ cTx

πTA ≤ cT & x ≥ 0

Minimize Maximize

▆

7

Solving Linear Programs

Three types of algorithms are available

Primal simplex algorithms (Dantzig 1947) Dual simplex algorithms (Lemke 1954)

Developed in context of game theory

Primal-dual log barrier algorithms

Interior-point algorithms (Karmarkar 1989) Reference: Primal-Dual Interior Point Methods, S.

Wright, 1997, SIAM

Primary focus: Dual simplex algorithms

8

Basic Solutions – Definition

Let B be an ordered set of m distinct indices (B1,…,Bm) taken from {1,…,n}. B is called a basis for (P) if AB is nonsingular. The variables xB are known as the basic variables and the variables xN as the non-basic variables, where N = {1,…,n}\B. The corresponding basic solution X∈ Rn is given by XN=0 and XB=AB-1 b. B is called (primal) feasible if XB ≥ 0.

Note: AX = b ⇒ ABXB + ANXN = b ⇒ AB XB = b ⇒ XB = AB

• 1b

9

Primal Simplex Algorithm

(Dantzig, 1947) Input: A feasible basis B and vectors XB = AB

• 1b and DN = cN – AN

TAB

• TcB.

Step 1: (Pricing) If DN ≥ 0, stop, B is optimal; else let

j = argmin{Dk : k∈N}.

Step 2: (FTRAN) Solve ABy=Aj. Step 3: (Ratio test) If y ≤ 0, stop, (P) is unbounded; else, let

i = argmin{XBk/yk: yk > 0}.

Step 4: (BTRAN) Solve AB

Tz = ei.

Step 5: (Update) Compute αN=-AN

• Tz. Let Bi=j. Update XB

(using y) and DN (using αN)

Note: xj is called the entering variable and xBi the leaving variable. The DN values are known as reduced costs – like partial derivatives

• f the objective function relative to the nonbasic variables.

10

Primal Simplex Example

11

The Simplex Algorithm

Consider the following simple LP:

Maximize 3x1 + 2x2 + 2x3 Subject to x1 + x3 ≤ 8 x1 + x2 ≤ 7 x1 + 2x2 ≤ 12 x1, x2, x3 ≥ 0

12

The Primal Simplex Algorithm

x1 x2 x3

(0,0,8) (0,6,8) (2,5,6) (2,5,0) (7,0,1) (7,0,0) (0,6,0)

Maximize z = 3x1 + 2x2 + 2x3

z = 0 z = 21 z = 23

Optimal

z = 28

Add slacks: Initial basis B = (4,5,6) Maximize 3x1 + 2x2 + 2x3 + 0x4 + 0x5 + 0x6 Subject to x1 + x3 + x4 = 8 x1 + x2 + x5 = 7 x1 + 2x2 + x6 = 12 x1, x2, x3,x4,x5,x6 ≥ 0

x1 enters, x5 leaves basis

D1 = rate of change of z relative to x1 = 21/7=3

13

Dual Simple Algorithm – Setup

Simplex algorithms apply to problems with constraints in equality form. We convert (D) to this form by adding the dual slacks d: Maximize bTπ Subject to ATπ + d = c π free, d ≥ 0 ⇔ ATπ ≤ c

14

Dual Simple Algorithm – Setup

Maximize bTπ Subject to ATπ + d = c π free, d ≥ 0 Given a basis B, the corresponding dual basic variables are π and dN. dB are the nonbasic

• variables. The corresponding dual basic

solution Π,D is determined as follows: DB=0 ⇒ Π = AB-TcB ⇒ DN=cN – ANTΠ B is dual feasible if DN ≥ 0.

= AB

T IB 0

AN

T 0 IN

π dB dN cB cN

15

Dual Simple Algorithm – Setup

Maximize bTπ Subject to ATπ + d = c π free, d ≥ 0

= AB

T IB 0

AN

T 0 IN

π dB dN cB cN

Observation: We may assume that every dual basis has the above form. Proof: Assuming that the primal has a basis is equivalent to assuming that rank(A)=m (# of rows), and this implies that all π variables can be assumed to be basic. This observation establishes a 1-1 correspondence between primal and dual bases.

▆

16

An Important Fact

If X and Π,D are corresponding primal and dual basic solutions determined by a basis B, then Π Tb = cTX. Hence, by weak duality, if B is both primal and dual feasible, then X is optimal for (P) and Π is

• ptimal for (D).

Proof: cTX = cBTXB (since XN=0) = Π TABXB (since Π = AB-TcB) = Π Tb (since ABXB=b)

▆

17

Dual Simplex Algorithm

(Lemke, 1954) Input: A dual feasible basis B and vectors XB = AB

• 1b and DN = cN – AN

TB-TcB.

Step 1: (Pricing) If XB ≥ 0, stop, B is optimal; else let

i = argmin{XBk : k∈{1,…,m}}.

Step 2: (BTRAN) Solve BTz = ei. Compute αN=-AN

Tz.

Step 3: (Ratio test) If αN ≤ 0, stop, (D) is unbounded; else, let

j = argmin{Dk/αk: αk > 0}.

Step 4: (FTRAN) Solve ABy = Aj. Step 5: (Update) Set Bi=j. Update XB (using y) and DN (using αN)

Note: dBi is the entering variable and dj is the leaving variable. (Expressed in terms of the primal: xBi is the leaving variable and xj is the entering variable)

18

Simplex Algorithms

Input: A primal feasible basis B and vectors XB=AB

• 1b & DN=cN – AN

TAB

• TcB.

Step 1: (Pricing) If DN ≥ 0, stop,

B is optimal; else, let j = argmin{Dk : k∈N}.

Step 2: (FTRAN) Solve ABy=Aj. Step 3: (Ratio test) If y ≤ 0, stop,

(P) is unbounded; else, let i = argmin{XBk/yk: yk > 0}.

Step 4: (BTRAN) Solve AB

Tz =

ei.

Step 5: (Update) Compute αN =

• AN
• Tz. Let Bi=j. Update XB

(using y) and DN (using αN)

Input: A dual feasible basis B and vectors XB=AB

• 1b & DN=cN – AN

TAB

• TcB.

Step 1: (Pricing) If XB ≥ 0, stop,

B is optimal; else, let i = argmin{XBk : k∈{1,…,m}}.

Step 2: (BTRAN) Solve AB

Tz =

• ei. Compute αN=-AN

Tz.

Step 3: (Ratio test) If αN ≤ 0,

stop, (D) is unbounded; else, let j = argmin{Dk/αk: αk > 0}.

Step 4: (FTRAN) Solve ABy =

Aj.

Step 5: (Update) Set Bi=j.

Update XB (using y) and DN (using αN)

19

Correctness: Dual Simplex Algorithm

Termination criteria

Optimality Unboundedness

Other issues

Finding starting dual feasible basis, or showing that no

feasible solution exists

Input conditions are preserved (i.e., that B is still a

feasible basis)

Finiteness

(DONE – by “An Important Fact” !!!)

20

Summary:

What we have done and what we have to do

Done

Defined primal and dual linear programs Proved the weak duality theorem Introduced the concept of a basis Stated primal and dual simplex algorithms

To do (for dual simplex algorithm)

Show correctness Describe key implementation ideas Motivation

21

Dual Unboundedness

(⇒ primal infeasible)

We carry out a key calculation As noted earlier, in an iteration of the dual The idea: Currently dBi = 0, and XBi < 0 has motivated

us to increase dBi to θ > 0, leaving the other components of dB at 0 (the object being to increase the

• bjective). Letting d,π be the corresponding dual

solution as a function of θ, we obtain dB = θ ei π = Π – θ z dN = DN – θ αN where αN and z are as computed in the algorithm.

dBi enters basis dj leaves basis in Maximize bTπ Subject to ATπ + d = c π free, d ≥ 0

22

(Dual Unboundedness – cont.)

Letting d,π be the corresponding dual solution as a

function of θ. Using αN and z from dual algorithm, dB = θ ei dN = DN – θ αN π = π – θ z.

Using θ > 0 and XBi < 0 yields

new_objective = πT b = (π – θ z)T b = πT b – θ XBi = old_objective – θ XBi > old_objective

Conclusion 1: If αN ≤ 0, then dN ≥ 0 ∀ θ > 0 ⇒ (D)

is unbounded.

Conclusion 2: If αN not≤ 0, then

dN ≥ 0 ⇒ θ ≤ Dj /αj ∀ αj > 0 ⇒ θmax = min{Dj /αj: αj > 0}

23

(Dual Unboundedness – cont.)

Finiteness: If DB > 0 for all dual feasible bases B, then

the dual simplex algorithm is finite: The dual

• bjective strictly increases at each iteration ⇒ no basis

repeats, and there are a finite number of bases.

There are various approaches to guaranteeing finiteness

in general:

Bland’s Rules: Purely combinatorial, bad in practice. CPLEX: A perturbation is introduced to guarantee

DB > 0.

24

Computational History

• f

Linear Programming

25

“A certain wide class of practical problems appears to be just

beyond the range of modern computing machinery. These problems occur in everyday life; they run the gamut from some very simple situations that confront an individual to those connected with the national economy as a whole. Typically, these problems involve a complex of different activities in which one wishes to know which activities to emphasize in order to carry out desired objectives under known limitations.”

George B. Dantzig, 1948

26

Application of LP & MIP - I

• Transportation-airlines
• Fleet assignment
• Crew scheduling
• Ground personnel scheduling
• Yield management
• Fuel allocation
• Passenger mix
• Booking control
• Maintenance scheduling
• Load balancing/freight packing
• Airport traffic planning
• Gate scheduling/assignment
• Upset recover and management
• Transportation-other
• Vehicle routing
• Freight vehicle scheduling and

assignment

• Depot/warehouse location
• Freight vehicle packing
• Public transportation system
• peration
• Rental car fleet management
• Process industries
• Plant production scheduling and

logistics

• Capacity expansion planning
• Pipeline transportation planning
• Gasoline and chemical blending

27

Application of LP & MIP - II

• Financial
• Portfolio selection and optimization
• Cash management
• Synthetic option development
• Lease analysis
• Capital budgeting and rationing
• Bank financial planning
• Accounting allocations
• Securities industry surveillance
• Audit staff planning
• Assets/liabilities management
• Unit costing
• Financial valuation
• Bank shift scheduling
• Consumer credit delinquency

management

• Check clearing systems
• Municipal bond bidding
• Stock exchange operations
• Debt financing
• Manufacturing
• Product mix planning
• Blending
• Manufacturing scheduling
• Inventory management
• Job scheduling
• Personnel scheduling
• Maintenance scheduling and planning
• Steel production scheduling
• Coal Industry
• Coal sourcing/transportation logistics
• Coal blending
• Mining operations management
• Forestry
• Forest land management
• Forest valuation models
• Planting and harvesting models

28

Application of LP & MIP - III

• Agriculture
• Production planning
• Farm land management
• Agricultural pricing models
• Crop and product mix decision models
• Product distribution
• Public utilities and natural resources
• Electric power distribution
• Power generator scheduling
• Power tariff rate determination
• Natural gas distribution planning
• Natural gas pipeline transportation
• Water resource management
• Alternative water supply evaluation
• Water reservoir management
• Public water transportation models
• Mining excavation models
• Oil and gas exploration and

production

• Oil and gas production scheduling
• Natural gas transportation scheduling
• Communications and computing
• Circuit board (VLSI) layout
• Logical circuit design
• Magnetic field design
• Complex computer graphics
• Curve fitting
• Virtual reality systems
• Computer system capacity planning
• Office automation
• Multiprocessor scheduling
• Telecommunications scheduling
• Telephone operator scheduling
• Telemarketing site selection

29

Application of LP & MIP - IV

• Food processing
• Food blending
• Recipe optimization
• Food transportation logistics
• Food manufacturing logistics and

scheduling

• Health care
• Hospital staff scheduling
• Hospital layout
• Health cost reimbursement
• Ambulance scheduling
• Radiation exposure models
• Pulp and paper industry
• Inventory planning
• Trim loss minimization
• Waste water recycling
• Transportation planning
• Textile industry
• Pattern layout and cutting optimization
• Production scheduling
• Government and military
• Post office scheduling and planning
• Military logistics
• Target assignment
• Missile detection
• Manpower deployment
• Miscellaneous applications
• Advertising mix/media scheduling
• Pollution control models
• Sales region definition
• Sales force deployment

30

CPLEX Across the World 1194 Cities – Excluding I SV Deployments

31

LP History

George Dantzig, 1947

Introduced LP and recognized it as more than a conceptual tool: Computing

answers important.

Invented “primal” simplex algorithm. First LP solved: Laderman, 9 cons., 77 vars., 120 MAN-DAYS.

First computer code – 1951 LP used commercially – Early 60s Powerful mainframe codes introduced – Early 70s Computational progress stagnated – Mid 80s Remarkable progress last 15 years (PCs, new computer

science and mathematics)

We now have three algorithms: Primal & Dual Simplex, Barrier

32

Example: A Production Planning Model

401,640 constraints 1,584,000 variables

Solution time line (2.0 GHz P4):

1988 (CPLEX 1.0): Houston, 13 Nov 2002

33

Example: A Production Planning Model

401,640 constraints 1,584,000 variables

Solution time line (2.0 GHz P4):

1988 (CPLEX 1.0):

8.0 days (Berlin, 21 Nov)

34

Example: A Production Planning Model

401,640 constraints 1,584,000 variables

Solution time line (2.0 GHz P4):

1988 (CPLEX 1.0):

15.0 days (Dagstuhl, 28 Nov)

35

Example: A Production Planning Model

401,640 constraints 1,584,000 variables

Solution time line (2.0 GHz P4):

1988 (CPLEX 1.0):

19.0 days (Amsterdam, 2 Dec)

36

Example: A Production Planning Model

401,640 constraints 1,584,000 variables

Solution time line (2.0 GHz P4):

1988 (CPLEX 1.0):

23.0 days (Houston, 6 Dec)

37

Example: A Production Planning Model

401,640 constraints 1,584,000 variables

Solution time line (2.0 GHz P4):

1988 (CPLEX 1.0):

29.8 days

Speedup: >43500x

1997 (CPLEX 5.0):

1.5 hours

2002 (CPLEX 8.0):

86.7 seconds

2003 (February):

59.1 seconds

38

BIG TEST: The testing methodology

Not possible for one test to cover 10+ years:

Combined several tests.

The biggest single test:

Assembled 680 real LPs Test runs: Using a time limit (4 days per LP) two chosen

methods would be compared as follows:

Run method 1: Generate 680 solve times Run method 2: Generate 680 solve times Compute 680 ratios and form GEOMETRIC MEAN (not

arithmetic mean!)

39

(part of a ) BIG Test Set

Model Rows Cols NZs Model Rows Cols NZs Model Rows Cols NZs

M1 16223 28568 88340 M28 31770 272372 829040 M55 99578 326504 2102273 M2 16768 39474 203112 M29 33440 56624 161831 M56 105127 154699 358171 M3 17681 165188 690273 M30 34994 87510 208179 M57 108393 112955 602948 M4 18262 23211 136324 M31 35519 43582 557466 M58 118158 487427 974854 M5 19103 33490 276895 M32 35645 34675 208769 M59 123964 93288 459680 M6 19374 180670 5392558 M33 36400 92878 246006 M60 125211 159109 457198 M7 19519 45832 124280 M34 38782 261079 1508199 M61 129181 467192 1025706 M8 19844 55528 152952 M35 39951 125000 381259 M62 155265 377918 930166 M9 19999 85191 170369 M36 41340 64162 370839 M63 175147 358239 1211488 M10 21019 115761 728432 M37 41344 163569 1928534 M64 179080 707556 1570514 M11 22513 99785 337746 M38 41366 78750 2110518 M65 185929 189867 2787708 M12 22797 63995 172018 M39 43387 107164 189864 M66 186441 23732 397080 M13 23610 44063 154822 M40 43687 164831 722066 M67 209760 363092 1061495 M14 23700 23005 169045 M41 44150 200077 4966017 M68 269640 1205640 6481640 M15 23712 31680 81245 M42 44211 37199 321663 M69 280756 920198 5936426 M16 24377 46592 2139096 M43 47423 81915 228565 M70 319256 638512 1231403 M17 26618 38904 1067713 M44 48548 163200 617683 M71 344297 559428 1909649 M18 27349 97710 288421 M45 54447 326504 1807146 M72 589250 1533590 5327318 M19 27441 15128 96118 M46 55020 117910 391081 M73 716772 1169910 2511088 M20 27899 26243 261968 M47 55463 191233 840986 M74 1000000 1685236 3370472 M21 28240 55200 161640 M48 60384 100078 485414 M75 1128152 1587664 4496138 M22 28420 164024 505253 M49 63856 144693 717229 M76 1285665 1313795 4998121 M23 29002 111722 2632880 M50 66185 157496 418321 M77 1372333 2627821 13321433 M24 29017 20074 2001102 M51 67745 111891 305125 M78 1709857 1903725 4959650 M25 29147 9984 1013168 M52 69418 612608 1722112 M79 5034171 7365337 25596099 M26 29724 98124 196524 M53 84840 316800 1899600 M80 5822606 15228380 30960543 M27 30190 57000 623730 M54 95011 197489 749771 M81 6662791 9781747 34053329 M82 10810259 19916187 49228262

40

LP Progress: 1988 – Present

(Operations Research, Jan 2002, pp. 3—15)

Algorithms (machine independent): Primal versus best of Primal/Dual/Barrier

3300x

Machines (workstations →PCs):

1600x

NET: Algorithm × Machine 5 300 000x (2 months/5300000 ~= 1 second)

41

Dual simplex vs. primal:

Dual 2.70x faster

Dual simplex vs. barrier:

Dual 1.06x faster

Algorithm comparison (≥ 50000 rows)

Algorithm comparison and other remarks …

(≥ 50000 rows)