CSDN博客

img yuanqingfei

线性规划的源程序(c语言版)

发表于2004/6/29 16:09:00  2377人阅读

分类: ☆Academic Endless☆

BBS 水木清华站 -- 精华区文章阅读


--------------------------------------------------------------------------------
发信人: FangQ (F the world), 信区: NumComp      
标  题: Re: 急寻线性规划的源程序(c语言版的)
发信站: BBS 水木清华站 (Sat Jun 30 15:45:48 2001)


Message 5 in thread
寄件者:Vic Smyth (vicsmyth@megsinet.net)
主旨:Re: Linear Programming
新闻群组:comp.programming
日期:1999/05/25

Tim,

    Don't know about C++, but Linear Programming and Simplex are covered
in chapter 10 of Numerical Recipies in C.  You can review the code
on-line at http://beta.ulib.org/webRoot/Books/Numerical_Recipes/

   -- Bill Clay


 
A program that we use at college in the linear programming (math) class is
called Lindo. But I don't know if it's available as freeware. Math software
Maple can also do Simplex. (Though I haven't taken the class yet and haven't
tried it.)


【 在 aliali (承轩) 的大作中提到: 】
: 谢谢。


--


        http://fangqq.dhs.org/
                    


※ 来源:·BBS 水木清华站 smth.org·[FROM: 129.170.67.237]
发信人: FangQ (F the world), 信区: NumComp      
标  题: Re: 急寻线性规划的源程序(c语言版的)
发信站: BBS 水木清华站 (Sat Jun 30 15:48:32 2001)

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.


Regards,

Stephen

=====================================================================
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<smallpos) then
                           begin
                              smallpos := value;
                              smallptr := basis[ row ];
                           end;
                        if (value<zero) and (value>largeneg) then
                           begin
                              largeneg := value;
                              largeptr := basis[ row ];
                           end;
                     end;
               end;
            if (row1<=nltcon+neqcon) then
               begin
                  write(row1:3);
                  if (smallpos <> mxval) then
                      write((getmatrix( row1, 0 )-smallpos):18:8)
                  else
                      write('No Limit':18);
                  write((getmatrix( row1, 0 )         ):18:8);
                  if (largeneg <> -mxval) then
                     write((getmatrix( row1, 0 )-largeneg):18:8)
                  else
                     write('No Limit':18);
                  if (smallptr <> 0) then
                     write(smallptr:6)
                  else
                     write('None':6);
                  if (largeptr <> 0) then
                     write(largeptr:6)
                  else
                     write('None':6);
               end
            else
               begin
                  write(row1:3);
                  if (largeneg <> -mxval) then
                     write((getmatrix( row1, 0 )+largeneg):18:8)
                  else
                     write('No Limit':18);
                  write((getmatrix( row1, 0 )         ):18:8);
                  if (smallpos <> mxval) then
                     write((getmatrix( row1, 0 )+smallpos):18:8)
                  else
                     write('No Limit':18);
                  if (largeptr <> 0) then
                     write(largeptr:6)
                  else
                     write('None':6);
                  if (smallptr <> 0) then
                     write(smallptr:6)
                  else
                     write('None':6);
               end;
            writeln;
         end;
      writeln(' ');
      writeln('* Sensitivity Analysis - Objective Coefficient Ranging *');
      writeln('Var':3,'Lower':18,'Current':18,'Upper':18,'InLo':6,'InUp':6);
      writeln('---':3,'-----':18,'-------':18,'-----':18,'----':6,'----':6);
      for col := 1 to nvar do
         begin
            smallpos := +mxval;
            smallptr := 0;
            largeneg := -mxval;
            largeptr := 0;
            if (basisp[ col ] = 0) then
             if (minmax < 0) then
               begin
                  smallpos := -minmax*getmatrix( trows2, col );
                  smallptr := col;
               end
              else
               begin
                  largeneg := -minmax*getmatrix( trows2, col );
                  largeptr := col;
               end
            else
             for row := 1 to tcols3 do
              if (basisp[ row ] = 0) then
               begin
                  value := getmatrix( basisp[ col ], row );
                  if ( value <> zero ) then
                     begin
                        value := minmax*getmatrix( trows2, row ) / value;
                        if (value>zero) and (value<smallpos) then
                           begin
                              smallpos := value;
                              smallptr := row;
                           end;
                        if (value<zero) and (value>largeneg) then
                           begin
                              largeneg := value;
                              largeptr := row;
                           end;
                     end;
               end;
            write(col:3);
            if (largeneg <> -mxval) then
               write((minmax*getmatrix( 0, col )+largeneg):18:8)
            else
               write('No Limit':18);
            write((minmax*getmatrix( 0, col )         ):18:8);
            if (smallpos <> mxval) then
               write((minmax*getmatrix( 0, col )+smallpos):18:8)
            else
               write('No Limit':18);
            if (largeptr <> 0) then
               write(largeptr:6)
            else
               write('None':6);
            if (smallptr <> 0) then
               write(smallptr:6)
            else
               write('None':6);
            writeln;
         end;
end;

procedure print;
var row, col : integer;
begin
   writeln;
   for row := 1 to trows2 do
      begin
         if (row > 0) and (row <= ncon) then
            write(basis[ row ]:2,' ')
         else
            write('   ');
         for col := 1 to tcols2 do
            write(getmatrix( row, col ):9:3,' ');
         writeln;
      end;
   writeln;
end;

begin
   initmatrix;
   writeln;
   writeln('*** Linear Programming - Simplex Algorithm ***');
   writeln;
   reader( error );
   problem;
   if ( error = 0 ) then
      simplex( error );
   if ( error = 0 ) or ( error = 1 ) then
      answer;
   if ( error < 0 ) then
      writeln('-- Inconsistent Data - Not Run --');
   if ( error = 2 ) then
      writeln('-- The Solution is Unbounded --');
   if ( error = 3 ) then
      writeln('-- The Problem is Infeasible --');
end.


********************
*** SAMPLE INPUT ***
********************

Photocopy - Sensitivity Analysis - Page 329
-1 5 4
2 0 3
    1.00     0.00     0.00     0.00      2.00
    0.00     1.00     0.00     0.00      1.50
  300.00  1000.00    10.00  1500.00   2500.00
   20.00   110.00    10.00     0.00    200.00
   15.00    10.00     8.00    48.00    100.00
   50.00    80.00    15.00   180.00

*********************
*** SAMPLE OUTPUT ***
*********************

*** Linear Programming - Simplex Algorithm ***

Filename: pc329.lp
 Problem: Photocopy - Sensitivity Analysis - Page 329

Minimize:
       50.00000000 * Var#  1 +
       80.00000000 * Var#  2 +
       15.00000000 * Var#  3 +
      180.00000000 * Var#  4

Subject to:
  Constraint #  1...        1.00000000 * Var#  1 +
                            0.00000000 * Var#  2 +
                            0.00000000 * Var#  3 +
                            0.00000000 * Var#  4 <=         2.00000000
  Constraint #  2...        0.00000000 * Var#  1 +
                            1.00000000 * Var#  2 +
                            0.00000000 * Var#  3 +
                            0.00000000 * Var#  4 <=         1.50000000
  Constraint #  3...      300.00000000 * Var#  1 +
                         1000.00000000 * Var#  2 +
                           10.00000000 * Var#  3 +
                         1500.00000000 * Var#  4 >=      2500.00000000
  Constraint #  4...       20.00000000 * Var#  1 +
                          110.00000000 * Var#  2 +
                           10.00000000 * Var#  3 +
                            0.00000000 * Var#  4 >=       200.00000000
  Constraint #  5...       15.00000000 * Var#  1 +
                           10.00000000 * Var#  2 +
                            8.00000000 * Var#  3 +
                           48.00000000 * Var#  4 >=       100.00000000

 * Your Variables : 1 through 4
 * Slack Variables : 5 through 6
 * Surplus Variables : 7 through 9
 * Artificial Variables : 10 through 12

* Value of Objective Function:       335.23437500

* Variable Analysis *
Var             Value      Reduced Cost
---             -----      ------------
  1        0.00000000       -4.29687500
  2        1.50000000        0.00000000
  3        6.90104167        0.00000000
  4        0.62065972        0.00000000

* Constraint Analysis *
Row         RHS Value     Slack/Surplus      Shadow Price
---         ---------     -------------      ------------
  1        2.00000000        2.00000000        0.00000000
  2        1.50000000        0.00000000       -0.46875000
  3     2500.00000000        0.00000000       -0.06250000
  4      200.00000000       34.01041667        0.00000000
  5      100.00000000        0.00000000       -1.79687500

* Sensitivity Analysis - Right Hand Side Ranging *
Row             Lower           Current             Upper OutLo OutUp
---             -----           -------             ----- ----- -----
  1        0.00000000        2.00000000          No Limit     5  None
  2        1.25469572        1.50000000        2.40506329     8     4
  3     1606.25000000     2500.00000000     3316.25000000     4     8
  4          No Limit      200.00000000      234.01041667  None     8
  5       73.88000000      100.00000000      815.00000001     8     4

* Sensitivity Analysis - Objective Coefficient Ranging *
Var             Lower           Current             Upper  InLo  InUp
---             -----           -------             -----  ----  ----
  1       45.70312500       50.00000000          No Limit     1  None
  2          No Limit       80.00000000       80.46875000  None     6
  3        1.20000000       15.00000000       15.16363636     9     6
  4      179.31645570      180.00000000      202.00000000     6     1

 

 


【 在 aliali (承轩) 的大作中提到: 】
: 谢谢。


--


        http://fangqq.dhs.org/
                    


※ 来源:·BBS 水木清华站 smth.org·[FROM: 129.170.67.237]
 

--------------------------------------------------------------------------------
[返回上一页]

0 0

相关博文

我的热门文章

img
取 消
img