Linear Programming in F#
Celso Axelrud
Jan/30/2011
Introduction
This document presents a Linear Programming function developed in F# and executed in the interactive environment.
I am targeting small LP programs on very fast computers.
I selected the simplest and shortest code from all references. It means that has no linear algebra operation and not subroutine calls.
The F# code presented here is based on the subroutine Linpro.
Fortran version at:
http://www.sie.arizona.edu/faculty/addenda/yak/SIE270/for/linpro
Pascal version at:
http://www.sie.arizona.edu/faculty/addenda/yak/SIE270/pas/linpro.pas
There are other references to a subroutine Linpro but it is not the same code.
Fortran version at:
http://www.slac.stanford.edu/accel/ilc/codes/lcopt/source/linpro.f
Python version at:
http://adorio-research.org/wordpress/?p=194
Source Code and Results
(*
Linpro - this function computes the solution of a linear programming
problem by using the Simplex method. The Objective function
is maximized. Parameters a,b, and c are presumed transmitted
as global variables declared in the calling program
Usage:
Linpro n m1 m2 m3
Parmeters:
a = matrix of coefficients (global)
b = right hand side vector (global)
c = coefficients of the objective function(global)
Input:
n = number of variables
m1 = number of constraints of the Le type
m2 = number of constraints of the Eq type
m3 = number of constraints of the Ge type
Output:
x.[i] = i-th component of the optimal solution
for i = 1, ... ,N
x.[n+1] = optimal value of the objective function
Kod = key showing the category of the Lp problem
if Kod = 1 then an optimal exists
if Kod = 2 then no feasible solution exists
if Kod = 3 then the objective function is not
bounded
Remarks:
matrix a needs to be re-assign before calling the function again
Conditions:
M1+M2+M3 must be less then or equal to 50
N must not be larger then 50 unless the subroutine is
redimenstion
Example:
Maximize Z = f(x,y) = 3x + 2y
subject to: 2x + y ≤ 18
2x + 3y ≤ 42
3x + y ≤ 24
x ≥ 0 , y ≥ 0
result:
(x,y) = (3,12)
Z=33
*)
let a:float[,]= Array2D.zeroCreate 52 151
let b:float[]= Array.zeroCreate 50
let c:float[]= Array.zeroCreate 50
//Linpro
let Linpro n m1 m2 m3 =
let x:float array = Array.zeroCreate 50
let mutable i,i0,i1,i2=0,0,0,0
let mutable j,j0,j1,j2,L=0,0,0,0,0
let mutable L1,L2,L3,Kod=0,0,0,0
let mutable m4,m5,mm,n1,n5,nm,nn=0,0,0,0,0,0,0
let mutable key=0
let mutable ifirst=0
let mutable u=0.0
let mutable Key=0
let eps=1.0e-3
let kk:int array = Array.zeroCreate 150
//Load up the initial tableau with slack and artificial
//variables, original and secondary objective functions
L1<- n + 1
L2<- m1+m2+m3+1
n1<- n + 1
//set the size of the appended matrix
mm<- m1+m2+m3+2
nn<- n+m1+m2+(2*m3)+1
//initialize both objectives to zero
for i= L2 to mm do
for j= 1 to nn do
a.[i,j]<- 0.0
L3<- m1+m2+m3
// initialize all slack and artificial coefficients to zero
for i= 1 to L3 do
for j= L1 to nn do
a.[i,j]<- 0.0
L<-n
if m1 > 0 then
// Fix slack coefficients for Le constraints
for i= 1 to m1 do
a.[i,L+i]<- 1.0
L<- L + m1
if m3 > 0 then
// Fix slack coefficients for Ge constraints
for i= 1 to m3 do
i1<-m1+m2+i
a.[i1,L+i]<- -1.0
L<- L + m3
m4<- m2 + m3
if (m4 > 0) then
// Fix artificial coefficient for for Eq and Ge constraints
for i= 1 to m4 do
a.[m1+i,L+i]<- 1.0
m5<- mm - 1
// Place the objective function
for j= 1 to n do
a.[m5,j]<- c.[j]
i2<- m1+m2+m3;
if (m4 > 0) then
i1<- m1 + 1
//construct the secondary objective function
for j= 1 to n do
for i= i1 to i2 do
a.[mm,j]<- a.[mm,j] + a.[i,j]
j1<- n+m1+1
j2<- n+m1+m3
if (m3 > 0) then
for j= j1 to j2 do
a.[mm,j]<- -1.0
for i= i1 to i2 do
a.[mm,nn]<- a.[mm,nn] + b.[i]
for i= 1 to i2 do
a.[i,nn]<- b.[i]
a.[m5,nn]<- 0.0
// Load up the initial basic solution
if (m1 > 0) then
for i= 1 to m1 do
kk.[i]<- i + n
if (m4 > 0) then
i1<- m1 + 1
for i= i1 to i2 do
kk.[i]<- n+m1+m3+i-i1+1
else
//If no constraint of Eq and Gt type, then no secondary
// objective function is needed
mm<- mm - 1
ifirst<- 0
u<- 0.0
while (ifirst = 0) || (u > eps) do
ifirst<- 1
n5<- nn - 1;
// Check for the pivot element find the largest
// coefficient of the objective
u<-a.[mm,1]
j0<- 1
for j= 1 to n5 do
if (a.[mm,j] > u) then
u<- a.[mm,j]
j0<- j
if (u > eps) then
Key<- 0
i0<- 0
while (Key = 0) && (i0 <= i2) do
i0<- i0 + 1
if (a.[i0,j0] >= eps) then Key<- 1
// Test if there is a positive a[i,j] in this column
// if not then the objective function is bounded
if (Key = 0) then
Kod<- 3
//goto 10;
else
if (i0 < i2) then
i1<- i0 + 1
u<- a.[i0,nn]/a.[i0,j0]
for i= i1 to i2 do
if (a.[i,j0] >= eps) then
if ((a.[i,nn]/a.[i,j0]) < u) then
i0<- i
u<- a.[i,nn]/a.[i,j0]
//Perform the elimination step
u<- a.[i0,j0]
for j= 1 to nn do
a.[i0,j]<- a.[i0,j]/u
for i= 1 to mm do
if (i <> i0) then
for j= 1 to nn do
if (j <> j0) then
a.[i,j]<- a.[i,j]-a.[i0,j]*a.[i,j0]
for i= 1 to mm do
if (i <> i0) then
a.[i,j0]<- 0.0;
//Register the newest basis vector and go back to
//perform the next step of elimination
kk.[i0]<- j0
else if (m4 > 0) then
L<- n+m1+m3;
// Check for feasible solution from the secondary
//objective
for i= 1 to i2 do
if (kk.[i] > L) then
Kod<- 2
//goto 10;
if Kod<>2 then
mm<- mm - 1
// Remove artificail variables and the secondary
//objective function
m4<- 0
nm<- nn-m2-m3
for i= 1 to mm do
a.[i,nm]<- a.[i,nn]
nn<- nm
ifirst<- 0
//From here we go back to continue elimination
if Kod<>2 || Kod<>3 then Kod<- 1
// Set up the optimal solution
for i= 1 to n do
x.[i]<- 0.0
for i= 1 to i2 do
j<- kk.[i];
x.[j]<- a.[i,nn]
x.[n+1]<- -a.[mm,nn]
x,Kod
(*
//Example 1
let n=2
let m1=3
let m2,m3=0,0
c.[1..2]<-[|3.0;2.0|]
b.[1..3]<-[|18.0;42.0;24.0|]
a.[1..3,1..2]<-array2D [[2.0;1.0];[2.0;3.0];[3.0;1.0]]
let r=Linpro n m1 m2 m3
printfn "LP category: %i" (snd r)
printfn "LP Objective Function: %A" ((fst r).[3])
printfn "LP Varaibles: %A" ((fst r).[1..2])
(*
>
LP category: 0
LP Objective Function: 33.0
LP Varaibles: [|3.0; 12.0|]
*)
*)
//Example 2
(*
R = –2x + 5y, subject to:
x <= 200
y <= 170
x >= 100
y >= 80
y + x >= 200
Solution: x=100,y=170
*)
let n=2
let m1=2
let m2=0
let m3=3
c.[1..2]<-[|-2.0;5.0|]
b.[1..3]<-[|200.0;170.0;100.0;80.0;200.0|]
a.Initialize()
a.[1..5,1..2]<-array2D [[1.0;0.0];[0.0;1.0];[1.0;0.0];[0.0;1.0];[1.0;1.0]]
let r=Linpro n m1 m2 m3
printfn "LP category: %i" (snd r)
printfn "LP Objective Function: %A" ((fst r).[3])
printfn "LP Varaibles: %A" ((fst r).[1..2])
