unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls, Grids, Buttons,math;
const long=250;
type
  TForm1 = class(TForm)
    panel2: TPanel;
    BitBtn2: TBitBtn;
    Panel1: TPanel;
    Label1: TLabel;
    Label2: TLabel;
    Edit1: TEdit;
    Edit2: TEdit;
    BitBtn1: TBitBtn;
    RadioGroup1: TRadioGroup;
    Label3: TLabel;
    StringGrid1: TStringGrid;
    BitBtn3: TBitBtn;
    procedure BitBtn1Click(Sender: TObject);
    procedure BitBtn2Click(Sender: TObject);
    procedure BitBtn3Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}
type  myset=set of 1..long;//自定義集合類型，集合元素類型為整型，
      //范圍1..long。
var m,n:integer; //m供應站的數量；n需求站的數量
    a:array[1..long] of Real;//供應站的供應量（以單位計算）
    b:array[1..long] of Real;//需求站的需要量（以單位計算）
    c:array[1..long] of array[1..long] of real;//運價矩陣
      //c[i,j]表示從第i個供應站到第j個需求站單位產品的運輸成本
    x:array[1..long] of array[1..long] of Real;//運輸矩陣
      //x[i,j]表示從第i個供應站供應給第j個需求站x[i,j]個單位的產品
    z:real;//目標函數值，即總的運輸成本
    r:Integer;//關鍵行，當前供應站的下標
    s:Integer;//關鍵列，當前需求站的下標
    rr:myset;//有待分配的供應站的下標的集合
    ss:myset;//有待分配的需求站的下標的集合

    u:array[1..long] of Real;//供應站的位勢值
    v:array[1..long] of Real;//需求站的位勢值
    d:array[1..long] of array[1..long] of real;//可能運費矩陣
    biaoshi:array[1..long] of array[1..long] of Char;//標識符
procedure chushihuabianliang;
var i,j:Integer;
begin
for i:=1 to long do
  begin
    for j:=1 to long do
      begin
        c[i,j]:=0;
        x[i,j]:=0;
      end;
    a[i]:=0;
    b[i]:=0;
  end;
rr:=[];
ss:=[];
z:=0;
end;
{-------------初始化變量-------------------------------------}

function  find_r:Integer;
var i:integer;
begin
Result:=-1;
for i:=1 to m do
  if i in rr then
     begin
       Result:=i;
       break;
     end;
end;
{--------------原理第一步，尋找待分配的供應站------------------------}

function  find_s:Integer;
var i:integer;
begin
Result:=-1;
for i:=1 to n do
  if i in ss then
     begin
       Result:=i;
       break;
     end;
end;
{--------------原理第一步，尋找待分配的需求站-------------------------}

procedure  first;
var i,j:Integer;
begin
for i:=1 to long do
  for j:=1 to long do
    d[i,j]:=-1;
for i:=1 to m do
  for j:=1 to n do
    if x[i,j]>0 then
      d[i,j]:=0;
end;
{--------------對應步近法求解的第一步-----------------------------------}

procedure  second;
var i,j:Integer;
    cishu:Integer;//累加while執行的次數
    temp_u:array[1..long] of Boolean;
    //長位long的boolean型變量數組，輔助判斷是否已求出所有的u
    temp_v:array[1..long] of Boolean;
    //長位long的boolean型變量數組，輔助判斷是否已求出所有的v
    temp_all:boolean;//判斷是否求出所有的u和v的位勢
    label 1;
begin
for i:=1 to long do
  begin
    temp_u[i]:=False;
    temp_v[i]:=False;
  end;//給局部變量賦初值為false表示沒有求出該點的位勢
for i:=1 to long do
  begin
    u[i]:=-1.11;
    v[i]:=-1.11;
  end;//初始化位勢值，此處無特別意義，只為防止浮點數運算的出錯
{------------初始化變量-----------------------------------------------}

1:v[1]:=0;
temp_v[1]:=True;
temp_all:=False;
cishu:=0;
while (not temp_all) and  (cishu<(m+n)) do
  begin
    for i:=1 to m do
      if not temp_u[i] then
      begin
        for j:=1 to n do
          if temp_v[j] and (d[i,j]=0) then
             begin
               u[i]:=c[i,j]-v[j];
               temp_u[i]:=True;
               break;
             end;
      end;//根據v的位勢,搜索求解u的位勢
    for j:=1 to n do
      if not temp_v[j] then
      begin
        for i:=1 to m do
          if temp_u[i] and (d[i,j]=0) then
            begin
               v[j]:=c[i,j]-u[i];
               temp_v[j]:=True;
               break;
             end;
      end;//根據u的位勢,搜索求解v的位勢
    temp_all:=True;
    for i:=1 to m do
      if temp_u[i]=False  then
         temp_all:=False;
    for j:=1 to n do
       if temp_v[j]=False  then
         temp_all:=False;
    cishu:=cishu+1;
  end;//while語句的結束

if cishu=2*m*n then
   begin
     //showmessage('退化解');
     for i:=1 to m do
       if temp_u[i]=False then
         for j:=1 to n do
           if temp_v[j]=False then
              begin
                d[i,j]:=0;
                goto  1;
              end;
   end;
end;
{--------------對應步近法求解的第二步-----------------------------------}

procedure  third;
var i,j:Integer;
begin
for i:=1 to m do
  for j:=1 to n do
    if d[i,j]=-1 then
       d[i,j]:=u[i]+v[j]-c[i,j];
  //計算非基變量的檢驗數
end;
{--------------對應步近法求解的第三步-----------------------------------}

procedure  four;
var i,j:Integer;
begin
for i:=1 to long do
  for j:=1 to long do
    if (x[i,j]=0) and (d[i,j]=0) then
        biaoshi[i,j]:='*';
end;
{--------------對應步近法求解的第四步-----------------------------------}

procedure five;
var i,j:Integer;
    temp:real; //最大值
begin
temp:=d[1,1]; //初始化最大值
r:=1;
s:=1;//記錄最大值的位置
for i:=1 to m do
  for j:=1 to n do
    if d[i,j]>temp then
      begin
        temp:=d[i,j];
        r:=i;
        s:=j;
      end;//尋找定位最大值
end;
{--------------對應步近法求解的第五步-----------------------------------}

procedure seven;
var i,j:Integer;
    temp:boolean;
    temp_s:integer;
    temp_r:integer;
    biaoshishuliang:Integer;
    label 1;
begin
biaoshi[r,s]:='+';//標明帶“+”的元素x[r,s]
temp_r:=r;//從第r行開始找
temp_s:=s;
biaoshishuliang:=1;
1:for j:=1 to n do
    if (j<>temp_s) and (x[temp_r,j]>0) then //找大于零的元素
       begin
         temp:=False;
         for i:=1 to m do
           if i<>temp_r then //除去當前行
           begin
             if  (x[i,j]>0) or (biaoshi[i,j]='*')  then
               begin
                 temp:=True;
                 Break;
               end;//判斷改列中是否有大于零或標識為*的元素
           end;
         if (biaoshishuliang>1) and odd(biaoshishuliang) and (j=s) then
            temp:=true;
         if temp then
           begin
              biaoshi[temp_r,j]:='-';
              temp_s:=j;
              biaoshishuliang:=biaoshishuliang+1;
              Break;
           end;
       end; //在r行中找出一個大于零的元素
temp:=False;
for i:=1 to m do
  if biaoshi[i,s]='-' then
    begin
      temp:=True;
      Break;
    end;
if not temp then
  begin
    for i:=1 to m do
      if (i<>temp_r) and ((x[i,temp_s]>0) or (biaoshi[i,temp_s]='*')) then
         begin
           temp:=False;
           for j:=1 to n do
             if  (j<>temp_s) and (x[i,j]>0) then
               begin
                 temp:=True;
                 Break;
               end;
           if temp then
             begin
               biaoshi[i,temp_s]:='+';
               temp_r:=i;
               biaoshishuliang:=biaoshishuliang+1;
               break;
             end;
         end;
     goto 1;
  end;
end;
{--------------對應步近法求解的第七步-----------------------------------}

procedure eight;
var i,j:Integer;
    temp:real; //最小值
    temp_s:integer;
    temp_r:integer;
begin
for i:=1 to m do
  for j:=1 to n do
    if biaoshi[i,j]='-' then
      begin
        temp:=abs(x[i,j]);
        temp_r:=i;
        temp_s:=j;
        break;
      end; //確定一個初始最小值
for i:=1 to m do
  for j:=1 to n do
    if (biaoshi[i,j]='-') and (abs(x[i,j])<temp) then
      begin
        temp:=abs(x[i,j]);
        temp_r:=i;
        temp_s:=j;
      end;//尋找最小值并定位
for i:=1 to m do
  for j:=1 to n do
     begin
       if  biaoshi[i,j]='+' then
          x[i,j]:=x[i,j]+temp
       else if biaoshi[i,j]='-' then
          x[i,j]:=x[i,j]-temp;
     end;//更新運輸矩陣
end;
{--------------對應步近法求解的第八步-----------------------------------}
procedure TForm1.BitBtn1Click(Sender: TObject);
var i,j:Integer;
begin
m:=StrToInt(Edit1.Text);
n:=StrToInt(Edit2.Text);
with StringGrid1 do
  begin
    RowCount:=m+2;
    ColCount:=n+2;
    for i:=0 to RowCount-1 do
      for j:=0 to ColCount-1 do
        Cells[j,i]:='';
    Width:=ColCount*(DefaultColWidth+2);
    Height:=RowCount*(DefaultRowHeight+2);
    top:=(panel2.Height-Height) div 2-20;
    Left:=(panel2.Width-Width) div 2;
    for i:=1 to m do
      Cells[0,i]:='A'+IntToStr(i);
    for j:=1 to n do
      Cells[j,0]:='B'+IntToStr(j);
    Cells[0,0]:='供\需';
    Cells[0,m+1]:='需求量';
    Cells[n+1,0]:='供應量';
  end;
BitBtn2.Enabled:=True;
StringGrid1.SetFocus;
end;

procedure TForm1.BitBtn2Click(Sender: TObject);
var i,j:Integer;
    temp:real;
    //臨時變量，輔助在供應站和需求站種尋找最小的供需量
begin
chushihuabianliang;
for i:=1 to m do
  for j:=1 to n do
    if StringGrid1.Cells[j,i]<>'' then
       c[i,j]:=StrToFloat(StringGrid1.Cells[j,i]);
//讀取運價矩陣
for i:=1 to m do
  if StringGrid1.Cells[n+1,i]<>'' then
     a[i]:=StrToFloat(StringGrid1.Cells[n+1,i]);
//讀取供應量
for j:=1 to n do
  if StringGrid1.Cells[j,m+1]<>'' then
     b[j]:=StrToFloat(StringGrid1.Cells[j,m+1]);
//讀取需求量
{----------讀取已知條件，運價矩陣、供應量、需求量------------}
for i:=1 to m  do
  rr:=rr+[i];//設置初始待分配的供應站的下標集合
for j:=1 to n do
  ss:=ss+[j];//設置初始待分配的需求站的下標集合

while (ss<>[]) or (rr<>[]) do
  begin
    r:=find_r;
    s:=find_s;//完成原理的第一步
    if r=-1 then
      begin
        temp:=b[s];
        b[s]:=b[s]-temp;
        if b[s]=0 then ss:=ss-[s];
      end
    else if s=-1 then
      begin
        temp:=a[r];
        a[r]:=a[r]-temp;
        if a[r]=0 then rr:=rr-[r];
      end
    else
    begin
    temp:=min(a[r],b[s]);
    x[r,s]:=temp;
    a[r]:=a[r]-temp;
    b[s]:=b[s]-temp;//完成原理的第二步
    if a[r]=0 then rr:=rr-[r];
    if b[s]=0 then ss:=ss-[s];
    end;
  end;
for i:=1 to m do
  for j:=1 to n do
    StringGrid1.Cells[j,i]:=FloatToStr(x[i,j]);
{----------求取并顯示最初可行的運輸方案-------------------------------}
BitBtn2.Enabled:=False;
BitBtn3.Enabled:=True;
end;

procedure TForm1.BitBtn3Click(Sender: TObject);
var i,j:Integer;
    label 1;
begin
1:first;
second;
third;
for i:=1 to m do
  for j:=1 to n do
    biaoshi[i,j]:=#0;
four;
five;
if d[r,s]<=0 then
  begin
    z:=0;
    for i:=1 to m do
      for j:=1 to n do
        begin
          StringGrid1.Cells[j,i]:=FloatToStr(x[i,j]);
          z:=z+x[i,j]*c[i,j];
        end;
     ShowMessage('最小運費是：'+FormatFloat('0.###',z));
     BitBtn3.Enabled:=False;
  end
else
  begin
    seven;
    eight;
    goto 1;
  end;
end;

end.