Here is a Pascal program to solve small problems using the simplex algorithm.
From: Stephen Gale
Newsgroups: sci.op-research
Subject: Re: SIMPLEX code for PC
Date: Sun, 5 May 1996 07:21:50 -0600
Attached below is a fairly simple Simplex Solver (written for
Turbo Pascal 3.0). Please let me know of how you use the
program, as well as any problems and any successes.
=====================================================================
Stephen F. Gale, B.Sc., CCP, I.S.P.
SFGale@FreeNet.Calgary.Ab.Ca
http://www.freenet.calgary.ab.ca/~sfgale/
=====================================================================
95/96 Program Director
Association for Systems Management - Calgary Chapter
http://www.freenet.calgary.ab.ca/populati/communit/asmcal/asmcal.html
=====================================================================
************************
*** SIMPLEX SOFTWARE ***
************************
program linearoptimization;
const rowmx = 72;
colmx = 112;
mxval = 1.0E+35;
zero = 0.0;
eqzero = 0.00000001;
var matrix : array [0..rowmx, 0..colmx] of real;
basis : array [1..rowmx] of integer;
basisp : array [1..rowmx] of integer;
minmax : real;
error : integer;
name : string[70];
filename : string[14];
ncon : integer; {number of constraints}
nvar : integer; {number of variables}
nltcon : integer; {number of less than constraints}
neqcon : integer; {number of equal to constraints}
ngtcon : integer; {number of greater than contraints}
trows1 : integer;
trows2 : integer;
tcols1 : integer;
tcols2 : integer;
tcols3 : integer;
function getmatrix( row, col : integer ) : real;
begin
getmatrix := matrix[ row, col ];
end;
procedure putmatrix( row, col : integer; value : real );
begin
matrix[ row, col ] := value;
end;
procedure initmatrix;
var row, col : integer;
begin
for row := 0 to rowmx do
for col := 0 to colmx do
putmatrix( row, col, zero );
end;
procedure price( var xcol : integer; trow : integer; var error : integer );
var quant, val : real;
col : integer;
begin
quant := -eqzero;
for col := 1 to tcols3 do
begin
val := getmatrix( trow, col );
if ( val < quant ) then
begin
xcol := col;
quant := val;
end;
end;
error := 0;
if ( quant = -eqzero ) then
error := 1;
end;
procedure leave( var xrow : integer; xcol : integer; var error : integer );
var quant, val : real;
row : integer;
begin
quant := mxval;
for row := 1 to ncon do
begin
val := getmatrix( row, xcol );
if ( val > eqzero ) then
begin
val := getmatrix( row, tcols2 ) / val;
if ( val < quant ) then
begin
xrow := row;
quant := val;
end;
end;
end;
error := 0;
if ( quant = mxval ) then
error := 2;
end;
procedure pivot( xrow, xcol : integer );
var value, val, vl : real;
row, col : integer;
begin
value := getmatrix( xrow, xcol );
for row := 1 to trows2 do
if ( row <> xrow ) then
begin
vl := getmatrix( row, xcol );
for col := 1 to tcols2 do
if ( col <> xcol ) then
begin
val := getmatrix( row, col ) - vl * getmatrix( xrow, col ) / value;
if ( abs( val ) < eqzero ) then
val := zero;
putmatrix( row, col, val );
end;
end;
for col := 1 to tcols2 do
putmatrix( xrow, col, getmatrix( xrow, col ) / value );
for row := 1 to trows2 do
putmatrix( row, xcol, zero );
putmatrix( xrow, xcol, 1.0);
basis[ xrow ] := xcol;
end;
procedure optimize( trow : integer; var error : integer);
var xrow, xcol, iterate : integer;
begin
repeat
price( xcol, trow, error );
if ( error = 0 ) then
leave( xrow, xcol, error );
if ( error = 0 ) then
pivot( xrow, xcol );
until ( error <> 0 )
end;
procedure simplex( var error : integer );
var val : real;
row, col : integer;
flag : boolean;
label 1000;
begin
if ( ncon <> nltcon ) then
begin
optimize( trows1, error );
if ( error > 1 ) then exit;
error := 3;
for row := 1 to ncon do
if ( basis[ row ] > tcols3 ) then
begin
if ( getmatrix( row, tcols2 ) > eqzero ) then exit;
flag := false;
col := 1;
repeat
if ( abs( getmatrix( row, col ) ) >= eqzero ) then
begin
pivot( row, col );
flag := true;
end;
col := col + 1;
until ( (flag) or (col > tcols3) );
end;
end;
error := 0;
optimize( trows2, error );
end;
procedure reader( var error : integer );
var row, col, column : integer;
value, amt : real;
filevar : text;
begin
error := 0;
writeln(con,'Problem file should be in the following format:');
writeln(con,' Line 1 : Up to 70 character problem description');
writeln(con,' Line 2 : +1 (for max) or -1 (for min); # of constraints; # of variables');
writeln(con,' Line 3 : # of <= constraints; # of = constraints; # of >= constraints');
writeln(con,' Next : Constraints coefficients and RHS value for each constraint');
writeln(con,' Last : Objective function coefficients');
writeln(con);
write (con,'Enter the filename containing the problem: ');
readln (con,filename);
assign(filevar,filename);
reset(filevar);
{ read the problem description }
readln( filevar, name );
{ read the minmax, number of constraints, number of variables }
readln( filevar, minmax, ncon, nvar );
minmax := -minmax;
{ read the number of less than, equal to, greater than contraints}
readln( filevar, nltcon, neqcon, ngtcon );
if ( ncon <> nltcon + neqcon + ngtcon ) then error := -1;
trows1 := ncon + 1;
trows2 := ncon + 2;
tcols1 := nvar + ncon + ngtcon;
tcols2 := tcols1 + 1;
tcols3 := nvar + nltcon + ngtcon;
{prepare matrix and basis}
for row := 1 to trows2 do
for col := 1 to tcols2 do
putmatrix( row, col, zero );
for row := 1 to ncon do
basis[ row ] := 0;
{prepare artificial and surplus variables}
for row := 1 to ncon do
if ( row <= nltcon ) then
begin
column := nvar + row;
basis[ row ] := column;
putmatrix( row, column, +1.0 );
end
else
begin
column := nvar + ngtcon + row;
basis[ row ] := column;
putmatrix( row, column, +1.0 );
if ( row > nltcon + neqcon ) then
begin
column := nvar - neqcon + row;
putmatrix( row, column, -1.0 );
putmatrix( trows1, column, +1.0 );
end
end;
{read matrix and right hand side}
for row := 1 to ncon do
begin
for col := 1 to nvar do
begin
read( filevar, value );
putmatrix( row, col, value );
end;
read( filevar, value );
putmatrix( row, 0, value );
putmatrix( row, tcols2, value );
readln( filevar );
end;
{ read the coefficients of the objective function }
for col := 1 to nvar do
begin
read( filevar, value);
putmatrix( 0, col, value * minmax );
putmatrix( trows2, col, value * minmax );
end;
readln( filevar );
{ calculate artifical variables }
for col := 1 to nvar do
begin
value := zero;
for row := nltcon+1 to ncon do
value := value - getmatrix( row, col );
putmatrix( trows1, col, value );
end;
close(filevar);
end;
procedure stats;
begin
writeln;
writeln(' * Your Variables : 1 through ', nvar);
if ( nltcon > 0 ) then
writeln(' * Slack Variables : ',nvar+1,' through ',nvar+nltcon);
if ( ngtcon > 0 ) then
writeln(' * Surplus Variables : ',nvar+nltcon+1,' through ',tcols3);
if ( nltcon <> ncon ) then
writeln(' * Artificial Variables : ',tcols3+1,' through ',tcols1);
end;
procedure setbasis;
var row, col : integer;
flag : boolean;
begin
for col := 1 to nvar+ncon do
begin
flag := false;
row := 1;
repeat
if ( basis[ row ] = col ) then
flag := true
else
row := row + 1;
until ( (flag) or (row > ncon) );
if (flag) then
basisp[ col ] := row
else
basisp[ col ] := 0;
end;
end;
procedure problem;
var row, col : integer;
begin
{filename and problem description}
writeln('Filename: ',filename);
writeln(' Problem: ',name);
writeln;
{objective function}
if minmax < 0 then
writeln('Maximize:')
else
writeln('Minimize:');
for col := 1 to nvar do
begin
write(minmax*getmatrix( trows2, col ):18:8,' * Var#',col:3);
if col <> nvar then
write(' + ');
writeln;
end;
writeln;
{constraints}
writeln('Subject to:');
for row := 1 to ncon do
begin
for col := 1 to nvar do
begin
if (col = 1) then
write(' Constraint #',row:3,'...')
else
write(' ':20);
write(getmatrix( row, col):18:8,' * Var#',col:3);
if col <> nvar then
writeln(' + ');
end;
if (row <= nltcon) then
write(' <= ')
else if (row > nltcon+neqcon) then
write(' >= ')
else
write(' = ');
writeln(getmatrix( row, tcols2 ):18:8);
end;
end;
procedure answer;
var smallpos : real;
smallptr : integer;
largeneg : real;
largeptr : integer;
value : real;
row, col, row1 : integer;
begin
setbasis;
stats;
{objective value}
writeln;
writeln('* Value of Objective Function: ',
-minmax*getmatrix( trows2, tcols2 ):18:8);
writeln;
writeln('* Variable Analysis *');
writeln('Var':3,'Value':18,'Reduced Cost':18);
writeln('---':3,'-----':18,'------------':18);
for col := 1 to nvar do
begin
write(col:3);
if (basisp[ col ] <> 0) then
write(getmatrix( basisp[ col ], tcols2 ):18:8)
else
write(0.0:18:8);
write(-minmax * getmatrix( trows2, col ):18:8);
writeln;
end;
writeln;
writeln('* Constraint Analysis *');
writeln('Row':3,'RHS Value':18,'Slack/Surplus':18,'Shadow Price':18);
writeln('---':3,'---------':18,'-------------':18,'------------':18);
for row := 1 to ncon do
begin
if (row <= nltcon) then
col := nvar + row
else if (row > nltcon+neqcon) then
col := nvar + row - neqcon
else
col := nvar + ngtcon + row;
write(row:3);
write(getmatrix( row, 0 ):18:8);
if (basisp[ col ] <> 0) then
write(getmatrix( basisp[ col ], tcols2 ):18:8)
else
write(0.0:18:8);
write(-minmax * getmatrix( trows2, col ):18:8);
writeln;
end;
writeln(' ');
writeln('* Sensitivity Analysis - Right Hand Side Ranging *');
writeln('Row':3,'Lower':18,'Current':18,'Upper':18,'OutLo':6,'OutUp':6);
writeln('---':3,'-----':18,'-------':18,'-----':18,'-----':6,'-----':6);
for row1 := 1 to ncon do
begin
if (row1 <= nltcon) then
col := nvar + row1
else if (row1 > nltcon+neqcon) then
col := nvar + row1 - neqcon
else
col := nvar + ngtcon + row1;
smallpos := +mxval;
smallptr := 0;
largeneg := -mxval;
largeptr := 0;
for row := 1 to ncon do
begin
value := getmatrix( row, col );
if ( value <> zero ) then
begin
value := getmatrix( row, tcols2 ) / value;
if (value>zero) and (value