Unregistered
|
Код писался 5 лет назад. Все забыто. Объявлений классов-матриц и векторов нет. Можно заменить на обычные массивы. Может разберешься. Пример внизу. | Код |
unit ччч;
interface
type
Float=double;
PFloat=^Float;
FloatArray=array[0..1000] of float;
PFloatArray=^FloatArray;
IntArray=array[0..1000] of integer;
PIntArray=^IntArray;
const
Big_IEEEDouble = 1.79E+0308; { max double }
Small_IEEEDouble = 4.94E-0324; { min double }
Tiny_IEEEDouble = 2.00E-0016; { min double such that 1+Tiny_IEEEDouble > 1 }
Big_IEEEExt = 5.95E+4931; { max extended }
Small_IEEEExt = 3.36E-4932; { min extended }
Tiny_IEEEExt = 2.00E-0019; { min extended such that 1+Tiny_IEEEextended > 1 }
function Simplex(A:TFloatMatrix; B:TFloatArray; C:TFloatArray; ConCnt,XCnt:integer;
EBasis:boolean; IBasis:TIntArray; X:TFloatArray; var OValue:extended; var Iter:longint;
MaxIter:longint; LECnt,GECnt:integer; Bi:TFloatMatrix; var Error:double):integer;
implementation
function Simplex(A:TFloatMatrix; B:TFloatArray; C:TFloatArray; ConCnt,XCnt:integer;
EBasis:boolean; IBasis:TIntArray; X:TFloatArray; var OValue:extended; var Iter:longint;
MaxIter:longint; LECnt,GECnt:integer; Bi:TFloatMatrix; var Error:double):integer;
{ Input:
A[ConCnt,XCnt] - coefficients of the constraints;
B[ConCnt] - right-hand side of the constraints;
C[XCnt] - vector of "costs" in the objective function;
ConCnt - dimention of array B (number of constraints);
XCnt - 2nd dimention of array A (number of variables);
EBasis - if false, the routine selects its own initial basis,
else indices of the initial basis should be in IBasis;
IBasis[ConCnt] - indices of the variables in the basis;
MaxIter - max number of iterations;
LECnt - number of <= constraints;
GECnt - number of >= constraints;
Output:
Result -
= 0 the problem was solved
= 1 the solution may have been obtained
= 2 MaxIter iterations were performed; more needed
= 3 sufficient accuracy could not be maintained to solve the problem
= 4 the problem has an unbounded solution
= 5 the problem has no solution
= 6 input error detected
X[] - dimention must be >=XCnt+LECnt+GECnt.
if Result in [0,1], it contains the values of the original,
slack and surplus variables;
OValue - if Result in [0,1], contains the value of the objective;
Iter - number of iterations used;
Bi[ConCnt,ConCnt] - the inverse of the basis matrix;
Error - the estimated relative error achieved.
}
label 22,30,32,100,200,300,360,430,440,450,460,500,580,600,601,610,630,650,660,
680,700,800;
var eps0,Errmn,Errmx,xmax,rtol,tmp,total,bnorm,binorm,t,dt,dsump,dsumn,s,tol:extended;
rmin,amax,rerr1,dsum,bmax,ratio,epsi,eps,cons,bmin,cmin,sgn,w:extended;
icount,mcheck,ms,m,i,j,ns,n,num,bflag,ibeg,iend,kj,l,ierr,iout,imin,lrow:integer;
ip,jp,nstep,nrow,ll,npos,ii,k:integer;
Xb,Y:PFloatArray; Basis,Indx:PIntArray;
procedure proc220;
var i,j:integer;
begin
for i:=0 to n-1 do X[i]:=0;
for i:=0 to m do begin
j:=IBasis[i];
if j<n then X[j]:=Xb^[i]
end
end;
procedure Done;
begin
FreeMem(Xb,ConCnt*SizeOf(float));
FreeMem(Y,ConCnt*SizeOf(float));
FreeMem(Basis,(ConCnt+Xcnt)*SizeOf(integer));
FreeMem(Indx,ConCnt*SizeOf(integer));
end;
procedure Croutl(Bi:TFloatArray; iend:integer; var ierr:integer);
label 20,30;
var i,maxdim,mcol,k,kcol,kk,nk,lk,l,lj,kj,nm1,km1,kp1,kcount,ik:integer;
kl,jcol,il,ncol,jj,j,nmk,lmin,n,m,ij,nmj:integer;
c,dsum,pmin,s:extended;
begin
n:=ConCnt;
m:=n-1;
maxdim:=sqr(n);
mcol:=iend*n;
if iend<>0 then begin
{ process the first iend columns of Bi }
kcol:=0;
for k:=0 to iend-1 do begin
kk:=kcol+k;
nk:=kcol+n-1;
lk:=kk;
while lk<=nk do begin
if Bi[lk]<0 then goto 20;
if Bi[lk]>0 then goto 30
end;
ierr:=1;
exit;
20:
lj:=mcol+lk-kcol;
while lj<maxdim do begin
Bi[lj]:=-Bi[lj];
inc(lj,n)
end;
30:
l:=lk-kcol;
Indx^[k]:=l;
if k<>l then begin
lj:=lk;
kj:=kk;
while kj<n do begin
{c:=Bi[kj]; Bi[kj]:=Bi[lj]; Bi[lj]:=c;}
Bi.SwapCells(kj,lj);
inc(lj,n);
inc(kj,n)
end
end;
inc(kcol,n)
end;
end;
{ process the remaining columns of Bi }
nm1:=n-2;
ierr:=0;
pmin:=0;
ibeg:=iend;
if ibeg<>n-1 then begin
k:=ibeg;
km1:=iend-1;
kp1:=k+1;
kcol:=mcol;
kk:=kcol+k;
for kcount:=ibeg to nm1 do begin
{ search for the k-th pivot element }
l:=k;
s:=abs(Bi[kk]);
for i:=kp1 to n-1 do begin
ik:=kcol+i;
c:=abs(Bi[ik]);
if s>=c then continue;
l:=i;
s:=c
end;
if (k>ibeg) and (s>=pmin) then
else begin
pmin:=s;
if s=0 then begin
ierr:=1;
exit
end
end;
{ interchanging rows k and l }
Indx^[k]:=l;
if k<>l then begin
kj:=mcol+k;
lj:=mcol+l;
while kj<maxdim do begin
{c:=Bi[kj]; Bi[kj]:=Bi[lj]; Bi[lj]:=c;}
Bi.SwapCells(kj,lj);
inc(lj,n);
inc(kj,n)
end
end;
{ compute the k-th row of U }
c:=Bi[kk];
if k<=ibeg then begin
kj:=kk+n;
while kj<maxdim do begin
Bi[kj]:=Bi[kj]/c;
inc(kj,n)
end
end
else begin
kl:=mcol+k;
for l:=ibeg to km1 do begin
Y^[l]:=Bi[kl];
inc(kl,n)
end;
kj:=kk+n;
while kj<maxdim do begin
jcol:=kj-k;
dsum:=-Bi[kj];
for l:=ibeg to km1 do begin
lj:=jcol+l;
dsum:=dsum+Y^[l]*Bi[lj]
end;
Bi[kj]:=-dsum/c;
inc(kj,n)
end;
end;
{ compute the k-th column of L }
km1:=k;
k:=kp1;
kp1:=k+1;
inc(kcol,n);
kk:=kcol+k;
for l:=ibeg to km1 do begin
lk:=kcol+l;
Y^[l]:=Bi[lk]
end;
for i:=k to n-1 do begin
il:=mcol+i;
dsum:=0;
for l:=ibeg to km1 do begin
dsum:=dsum+Bi[il]*Y^[l];
inc(il,n)
end;
Bi[il]:=Bi[il]-dsum
end
end;
end;
{ check the k-th pivot element }
ncol:=maxdim-n-1;
c:=abs(Bi[maxdim-1]);
if c>pmin then
else
if c=0 then begin
ierr:=1; exit
end;
{ replace L with the inverse of L }
if ibeg<>n-1 then begin
jj:=mcol+ibeg;
i:=n+1;
for j:=ibeg to nm1 do begin
Bi[jj]:=1/Bi[jj];
Y^[j]:=Bi[jj];
kj:=jj;
for km1:=j to nm1 do begin
k:=km1+1;
inc(kj);
dsum:=0;
kl:=kj;
for l:=j to km1 do begin
dsum:=dsum+Bi[kl]*Y^[l];
inc(kl,n)
end;
Bi[kj]:=-dsum/Bi[kl];
Y^[k]:=Bi[kj]
end;
inc(jj,i)
end
end;
Bi[maxdim-1]:=1/Bi[maxdim-1];
if n=1 then exit;
{ solve UX=Y where Y is the inverse of L }
for nmk:=0 to nm1 do begin
k:=n-nmk-1;
lmin:=MaxInt(ibeg,k);
kl:=lmin*n+k-1;
for l:=lmin to n-1 do begin
Y^[l]:=Bi[kl];
Bi[kl]:=0;
inc(kl,n)
end;
kj:=mcol+k-1;
while kj<maxdim do begin
dsum:=-Bi[kj];
lj:=(kj-k)+lmin+1;
for l:=lmin to n-1 do begin
dsum:=dsum+Y^[l]*Bi[lj];
inc(lj)
end;
Bi[kj]:=-dsum;
inc(kj,n)
end
end;
jcol:=ncol-n;
for nmj:=0 to nm1 do begin
j:=n-nmj-2;
k:=Indx^[j];
if j<>k then begin
ij:=jcol;
ik:=k*n;
for i:=0 to n-1 do begin
inc(ij);
{c:=Bi[ij]; Bi[ij]:=Bi[ik]; Bi[ik]:=c;}
Bi.SwapCells(ij,ik);
inc(ik)
end
end;
dec(jcol,n)
end;
end;
begin
eps0:=Tiny_IEEEExt;
Errmn:=10*eps0; { the smallest relative error tolerance used }
Errmx:=1E-5; { the largest relative error tolerance used }
xmax:=Big_IEEEExt; { the largest positive floating point number }
m:=ConCnt-1;
Iter:=0;
icount:=0;
mcheck:=MinInt(5,1+ConCnt div 15);
OValue:=0;
{ check for input errors }
ms:=LECnt+GECnt;
if (ConCnt<2) or (XCnt<2) or (ms>ConCnt) then begin
Result:=6; exit
end;
for i:=0 to m do
if B[i]<0 then begin
Result:=6; exit
end;
GetMem(Xb,ConCnt*SizeOf(float));
FillChar(Xb^,ConCnt*SizeOf(float),0);
GetMem(Y,ConCnt*SizeOf(float));
GetMem(Basis,(ConCnt+XCnt)*SizeOf(integer));
GetMem(Indx,ConCnt*SizeOf(integer));
try
rtol:=xmax;
for i:=0 to XCnt-1 do
if C[i]<>0 then rtol:=MinFloat(abs(C[i]),rtol);
rtol:=Errmx*rtol;
{ formation of the IBasis and Basis arrays }
ns:=XCnt+LECnt;
n:=ns+GECnt;
if not EBasis then goto 30;
num:=n;
for i:=0 to m do
if IBasis[i]>=n then inc(num);
goto 32;
22:
if not EBasis then begin
Result:=3;
Proc220;
exit
end;
EBasis:=False;
Result:=0;
30:
num:=XCnt+ConCnt;
for i:=0 to m do IBasis[i]:=XCnt+i;
32:
bflag:=0;
FillChar(Basis^[0],n*SizeOf(integer),0);
for i:=0 to m do Basis^[IBasis[i]]:=1;
if not EBasis then begin
{ calculation of Xb and Bi when EBasis=False }
Error:=Errmn;
Bi.Clear;
for j:=0 to m do begin
Xb^[j]:=B[j];
Bi[j,j]:=1
end;
if GECnt=0 then goto 630;
for j:=LECnt to ms-1 do begin
Xb^[j]:=-Xb^[j];
Bi[j,j]:=-1
end;
goto 601
end;
100:
{ reorder the basis }
ibeg:=0;
iend:=m;
for i:=0 to m do
if IBasis[i]>=XCnt then begin
Indx^[ibeg]:=IBasis[i];
inc(ibeg)
end
else begin
Indx^[iend]:=IBasis[i];
dec(iend)
end;
if iend=m then goto 22;
for i:=0 to m do IBasis[i]:=Indx^[i]; { IBasis[1:ConCnt]=Indx[1:ConCnt] }
{ reinversion of the basis matrix }
for j:=0 to m do begin
kj:=IBasis[j];
if kj<XCnt then begin
for i:=0 to m do Bi[i,j]:=A[i,kj]; { Bi[1:ConCnt,j]=A[1:ConCnt,kj] }
continue
end;
l:=kj-XCnt;
Bi.FillCol(0,m,j,0); { Bi[1:ConCnt,j]=0 }
if (kj<ns) or (kj>=n) then Bi[l,j]:=1 else Bi[l,j]:=-1
end;
icount:=0;
Croutl(TFloatArray(Bi),iend+1,ierr);
if ierr<>0 then goto 580;
{ check the accuracy of Bi and reset rerr }
bnorm:=0;
for j:=0 to m do begin
kj:=IBasis[j];
if kj<XCnt then begin
total:=0;
for i:=0 to m do total:=total+abs(A[i,kj])
end
else total:=1;
bnorm:=MaxFloat(bnorm,total);
end;
binorm:=0;
for j:=0 to m do begin
total:=0;
for i:=0 to m do total:=total+abs(Bi[i,j]);
binorm:=MaxFloat(binorm,total)
end;
Error:=MaxFloat(Errmn,eps0*bnorm*binorm);
if Error>1E-2 then goto 580;
bflag:=0;
{ recalculation of Xb }
for i:=0 to m do begin
dsump:=0; dsumn:=0;
for l:=0 to m do begin
dt:=Bi[i,l]*B[l];
if dt>0 then dsump:=dsump+dt else dsumn:=dsumn+dt;
end;
Xb^[i]:=dsump+dsumn;
s:=dsump;
t:=dsumn;
tol:=Errmx*MaxFloat(s,-t);
if abs(Xb^[i])<tol then Xb^[i]:=0
end;
goto 601;
200:
{ find the next vector A(--,jp) to be inserted into the basis }
jp:=-1;
rmin:=0;
if Nstep=3 then rmin:=-rtol;
for j:=0 to XCnt-1 do begin
if (Basis^[j]<>0) or (X[j]>=rmin) then continue;
jp:=j;
rmin:=X[j]
end;
if n<>XCnt then begin
rmin:=rmin*1.1;
for j:=XCnt to n-1 do begin
if (Basis^[j]<>0) or (X[j]>=rmin) then continue;
jp:=j;
rmin:=X[j]
end;
end;
if jp<>-1 then goto 300;
if nstep<2 then goto 800;
if nstep=2 then begin
{ completion of the nstep=2 case }
for i:=0 to m do begin
if IBasis[i]<n then continue;
if Xb^[i]>0 then goto 800
end;
goto 680;
end;
{ completion of the nstep=3 case }
if Error>1E-2 then begin
if icount>=5 then goto 100;
Result:=5{1};
end
else Result:=0;
goto 800;
{ if MaxIter iterations have not been performed then begin the next iteration.
compute the jp-th column of Bi*A and store it in Y }
300:
if iter>=MaxIter then begin
Result:=2;
Proc220;
exit
end;
inc(iter);
inc(icount);
if jp>=ns then begin
l:=jp-XCnt;
for i:=0 to m do Y^[i]:=-Bi[i,l]
end
else
if jp>=XCnt then begin
l:=jp-XCnt;
for i:=0 to m do Y^[i]:=Bi[i,l]
end
else begin
nrow:=-1;
amax:=0;
for i:=0 to m do begin
if A[i,jp]=0 then continue;
inc(nrow);
Indx^[nrow]:=i;
amax:=MaxFloat(abs(A[i,jp]),amax)
end;
if nrow=-1 then begin
Result:=4;
Proc220;
exit
end;
rerr1:=Errmx*amax;
for i:=0 to m do begin
dsum:=0;
for ll:=0 to nrow do begin
l:=Indx^[ll];
dsum:=dsum+Bi[i,l]*A[l,jp];
end;
Y^[i]:=dsum;
if abs(Y^[i])>=5E-3 then continue;
bmax:=0;
for l:=0 to m do bmax:=MaxFloat(abs(Bi[i,l]),bmax);
tol:=rerr1*bmax;
if abs(Y^[i])<tol then Y^[i]:=0
end;
end;
for i:=0 to m do
if Y^[i]<>0 then goto 360;
X[jp]:=0;
dec(Iter);
Dec(icount);
goto 200;
360:
if nstep=2 then goto 430;
if nstep>2 then goto 440;
{ finding the variable Xb[ip] to be made nonbasic for the nstep=1 case }
npos:=-1;
ip:=-1;
eps:=0;
epsi:=xmax;
for i:=0 to m do begin
if (Xb^[i]<0) or (Y^[i]<=0) then continue;
ratio:=Xb^[i]/Y^[i];
if ratio<epsi then begin
epsi:=ratio;
npos:=0;
Indx^[0]:=i;
continue
end
else
if ratio>epsi then continue;
inc(npos);
Indx^[npos]:=i
end;
if npos=-1 then begin
for i:=0 to m do begin
if (Xb^[i]>=0) or (Y^[i]>=0) then continue;
ratio:=Xb^[i]/Y^[i];
if ratio<eps then continue;
eps:=ratio;
ip:=i
end;
goto 500
end;
if epsi<>0 then begin
for i:=0 to m do begin
if (Xb^[i]>=0) or (Y^[i]>=0) then continue;
ratio:=Xb^[i]/Y^[i];
if (ratio>epsi) or (ratio<eps) then continue;
eps:=ratio;
ip:=i
end;
if ip<>-1 then goto 500
end;
460:
{ tie breaking procedure }
ip:=Indx^[0];
if npos=0 then goto 500;
ip:=-1;
bmin:=xmax; cmin:=xmax;
for ii:=0 to npos do begin
i:=Indx^[ii];
l:=IBasis[i];
if l<XCnt then begin
if C[l]<=0 then cmin:=MinFloat(0,cmin);
if C[l]>cmin then continue;
imin:=i;
cmin:=C[l]
end
else begin
if l>=n then begin
ip:=i;
goto 500
end;
lrow:=l-XCnt;
s:=B[lrow];
if lrow>=LECnt then begin
s:=-s;
bmin:=MinFloat(0,bmin);
end;
if s>bmin then continue;
ip:=i;
bmin:=s
end
end;
if (cmin<=0) or (ip=-1) then ip:=imin;
500:
{ transformation of Xb }
if Xb^[ip]<>0 then begin
cons:=Xb^[ip]/Y^[ip];
for i:=0 to m do begin
s:=Xb^[i];
Xb^[i]:=Xb^[i]-cons*Y^[i];
if Xb^[i]>=0 then continue;
if (s>=0) or (Xb^[i]>=Errmx*s) then Xb^[i]:=0
end;
Xb^[ip]:=cons
end;
{ transformation of Bi }
for j:=0 to m do begin
if Bi[ip,j]=0 then continue;
cons:=Bi[ip,j]/Y^[ip];
for i:=0 to m do Bi[i,j]:=Bi[i,j]-cons*Y^[i];
Bi[ip,j]:=cons
end;
{ updating IBasis and Basis }
iout:=IBasis[ip];
IBasis[ip]:=jp;
Basis^[iout]:=0;
Basis^[jp]:=1;
if iout>=n then dec(num);
{ check the accuracy of Bi and reset rerr }
if Error<=1E-2 then begin
k:=0;
for j:=0 to m do begin
kj:=IBasis[j];
if kj>=XCnt then continue;
total:=0;
for i:=0 to m do total:=total+Bi[j,i]*A[i,kj]; {Dot_Product(Bi[j,1:ConCnt],A[1:ConCnt,kj]);}
Error:=MaxFloat(Error,abs(1-total));
inc(k);
if k>=mcheck then break
end;
if Error<=1E-2 then goto 600;
end;
{ the accuracy criteria are not satisfied }
if icount<5 then goto 600;
bflag:=1;
goto 100;
580:
if iter=0 then begin
Result:=6; exit
end;
if bflag=0 then begin
Result:=3;
Proc220;
exit
end;
bflag:=0;
ip:=0;
while ip<=m do begin
if jp=IBasis[ip] then break;
inc(ip)
end;
IBasis[ip]:=iout;
Basis^[jp]:=0;
Basis^[iout]:=1;
if iout>=n then inc(num);
goto 100;
430:
{ finding the variable Xb[ip] to be made nonbasic for the nstep=2 case }
npos:=-1;
epsi:=xmax;
for i:=0 to m do begin
if Y^[i]<=0 then continue;
ratio:=Xb^[i]/Y^[i];
if ratio<epsi then begin
epsi:=ratio;
npos:=0;
Indx^[0]:=i
end
else
if ratio>epsi then continue;
inc(npos);
Indx^[npos]:=i
end;
goto 450;
440:
{ finding the variable Xb[ip] to be made nonbasic for the nstep=3 case }
npos:=-1;
epsi:=xmax;
for i:=0 to m do begin
if Y^[i]<0 then begin
if IBasis[i]<n then continue;
ip:=i;
goto 500
end
else
if Y^[i]>0 then begin
ratio:=Xb^[i]/Y^[i];
if ratio<epsi then begin
epsi:=ratio;
npos:=0;
Indx^[0]:=i
end
else
if ratio>epsi then continue;
inc(npos);
Indx^[npos]:=i
end;
end;
450:
if npos<>-1 then goto 460;
if icount>=5 then goto 100;
Result:=4;
Proc220;
exit;
600:
{ set up the R array for the nstep=1 case }
if nstep=2 then goto 630;
if nstep>2 then goto 700;
601:
for j:=0 to m do
if Xb^[j]<0 then goto 610;
goto 630;
610:
nstep:=1;
i:=-1;
for l:=0 to m do begin
if Xb^[l]>=0 then continue;
inc(i);
Indx^[i]:=l;
end;
for j:=0 to m do begin
dsump:=0; dsumn:=0;
for ll:=0 to i do begin
l:=Indx^[ll];
if Bi[l,j]<0 then dsumn:=dsumn+Bi[l,j]
else
if Bi[l,j]>0 then dsump:=dsump+Bi[l,j];
end;
Y^[j]:=dsump+dsumn;
s:=dsump;
t:=dsumn;
tol:=Errmx*MaxFloat(s,-t);
if abs(Y^[j])<tol then Y^[j]:=0
end;
goto 650;
630:
{ setup the R array for the nstep=2 case }
if n=num then goto 680;
nstep:=2;
i:=-1;
for l:=0 to m do begin
if IBasis[l]<n then continue;
inc(i);
Indx^[i]:=l
end;
for j:=0 to m do begin
dsump:=0; dsumn:=0;
for ll:=0 to i do begin
l:=Indx^[ll];
if Bi[l,j]<0 then dsumn:=dsumn+Bi[l,j]
else
if Bi[l,j]>0 then dsump:=dsump+Bi[l,j];
end;
Y^[j]:=-(dsump+dsumn);
s:=dsump;
t:=dsumn;
tol:=Errmx*MaxFloat(s,-t);
if abs(Y^[j])<tol then Y^[j]:=0
end;
650:
for j:=0 to XCnt-1 do
if Basis^[j]=0 then begin
tmp:=0;
for i:=0 to m do tmp:=tmp+Y^[i]*A[i,j];
X[j]:=tmp
end
else X[j]:=0;
660:
if XCnt<>ns then begin
for j:=XCnt to ns-1 do begin
X[j]:=0;
if Basis^[j]<>0 then continue;
X[j]:=Y^[j-XCnt]
end
end;
if ns=n then goto 200;
for j:=ns to n-1 do begin
X[j]:=0;
if Basis^[j]<>0 then continue;
X[j]:=-Y[j-XCnt]
end;
goto 200;
680:
{ setup the R array for the nstep=3 case }
nstep:=3;
for j:=0 to m do begin
dsum:=0;
for l:=0 to m do begin
ll:=IBasis[l];
if ll<XCnt then dsum:=dsum+C[ll]*Bi[l,j]
end;
Y^[j]:=dsum
end;
for j:=0 to XCnt-1 do begin
X[j]:=0;
if Basis^[j]<>0 then continue;
dsum:=-C[j];
for i:=0 to m do dsum:=dsum+Y^[i]*A[i,j];
X[j]:=dsum;
if dsum>=0 then continue;
tol:=Errmx*abs(C[j]);
if abs(X[j])<tol then X[j]:=0
end;
goto 660;
700:
{ update the R array for the nstep=3 case }
cons:=X[jp];
for j:=0 to XCnt-1 do begin
if Basis^[j]<>0 then X[j]:=0
else begin
total:=0;
for i:=0 to m do total:=total+Bi[ip,i]*A[i,j];
X[j]:=X[j]-cons*total;
if X[j]>=0 then continue;
tol:=Errmx*abs(C[j]);
if abs(X[j])<tol then X[j]:=0
end
end;
if XCnt<>ns then begin
for j:=XCnt to ns-1 do begin
if Basis^[j]<>0 then X[j]:=0 else X[j]:=X[j]-cons*Bi[ip,j-XCnt]
end;
end;
if ns=n then goto 200;
for j:=ns to n-1 do begin
if Basis^[j]<>0 then X[j]:=0 else X[j]:=X[j]+cons*Bi[ip,j-XCnt]
end;
goto 200;
{ refine Xb and stote the result in Y }
800:
FillChar(Y^,ConCnt*SizeOf(float),0);
ll:=-1;
for j:=0 to m do begin
kj:=IBasis[j];
if kj<XCnt then begin
inc(ll);
Indx^[ll]:=j;
continue
end;
l:=kj-XCnt;
if (kj<ns) or (kj>=n) then Y^[l]:=Xb^[j] else Y^[l]:=-Xb^[j]
end;
if ll=-1 then
for i:=0 to m do X[i]:=B[i]-Y^[i]
else
for i:=0 to m do begin
dsum:=Y^[i];
for ii:=0 to ll do begin
j:=Indx^[ii];
kj:=IBasis[j];
dsum:=dsum+A[i,kj]*Xb^[j]
end;
X[i]:=B[i]-dsum
end;
rerr1:=MinFloat(Errmx,Error);
for i:=0 to m do begin
Y^[i]:=0;
if Xb^[i]<0 then begin
sgn:=-1;
dsump:=0;
dsumn:=Xb^[i]
end
else
if Xb^[i]>0 then begin
sgn:=1;
dsump:=Xb^[i];
dsumn:=0
end
else continue;
for l:=0 to m do begin
dt:=Bi[i,l]*X[l];
if dt>0 then dsump:=dsump+dt else dsumn:=dsumn+dt;
end;
w:=dsump+dsumn;
if w=0 then continue;
if sgn<>ForSign(1,w) then continue;
s:=dsump;
t:=dsumn;
tol:=rerr1*MaxFloat(s,-t);
if abs(w)>tol then Y^[i]:=w
end;
case nstep of
1:begin
{ check the refinement (nstep=1) }
for i:=0 to m do begin
if Y^[i]<0 then begin
if Y^[i]<-Errmx then begin
if icount>=5 then goto 100;
Result:=5{1};
Proc220;
exit
end;
Y^[i]:=0
end;
Xb^[i]:=Y^[i]
end;
goto 630;
end;
2: begin
{ check the refinement (nstep=2) }
for i:=0 to m do begin
if IBasis[i]>=n then begin
if Y^[i]>Errmx then begin
if icount>=5 then goto 100;
Result:=5{1};
Proc220;
exit
end;
Y^[i]:=0
end;
Xb^[i]:=Y^[i]
end;
goto 680;
end;
else begin
{ compute OValue (nstep=3) }
dsum:=0;
for i:=0 to m do begin
j:=IBasis[i];
if j<XCnt then dsum:=dsum+C[j]*Y^[i];
Xb^[i]:=Y^[i]
end;
OValue:=dsum;
Proc220;
exit
end;
end; { case }
finally
Done
end
end;
end.
procedure TForm1.TestLinProg;
var i,j,Res,XDim:integer; ConCnt,XCnt,LECnt,GECnt:integer; EBasis:boolean;
A:TGFloatMatrix; B:TFloatArray; C:TFloatArray; IBasis:TIntArray; X:TFloatArray;
v,OValue:extended; Iter:longint; Error:double;
MaxIter:longint; Bi:TGFloatMatrix; s:string;
begin
(*
ConCnt:=2;
XCnt:=2;
LECnt:=2;
GECnt:=0;
*)
ConCnt:=SpinEdit2.Value;
XCnt:=SpinEdit1.Value;
LECnt:=SpinEdit3.Value;
GECnt:=SpinEdit4.Value;
EBasis:=False;
MaxIter:=100;
A:=TGFloatMatrix.Create(ConCnt,XCnt);
B:=TFloatArray.Create(ConCnt);
C:=TFloatArray.Create(XCnt);
IBasis:=TIntArray.Create(ConCnt);
XDim:=XCnt+LECnt+GECnt;
X:=TFloatArray.Create(XDim);
Bi:=TGFloatMatrix.Create(ConCnt,ConCnt);
for i:=0 to ConCnt-1 do begin
s:=Grid1.Cells[0,i+1];
val(s,v,Res);
B[i]:=v;
if i<XCnt then begin
s:=Grid2.Cells[0,i];
val(s,v,Res);
C[i]:=v;
end;
for j:=0 to XCnt-1 do begin
s:=Grid1.Cells[j+1,i+1];
val(s,v,Res);
A[i,j]:=v
end
end;
(*
A.InitFromConstArray([2,3,4,1{1,-1,1,0,1,4,0,1}{3,2,4,5}]);
B.InitFromConstArray([9,6{20,20,10,5}{1700,1600}]);
C.InitFromConstArray([6,2{3,4}{2,4}]);
*)
try
Res:=Simplex(A,B,C,ConCnt,XCnt,EBasis,IBasis,X,OValue,Iter,MaxIter,LECnt,GECnt,Bi,Error);
Memo1.Text:=IntToStr(Res);
Memo1.Lines.Add(FloatToStr(OValue));
Memo1.Lines.Add('Iters='+IntToStr(Iter));
for i:=0 to XCnt-1 do Memo1.Lines.Add('x['+IntToStr(i+1)+']='+FloatToStr(X[i]));
Memo1.Lines.Add('constrains:');
for i:=0 to ConCnt-1 do begin
v:=0;
for j:=0 to XCnt-1 do v:=v+A[i,j]*X[j];
Memo1.Lines.Add(FloatToStr(v)+' ('+FloatToStr(B[i])+')');
end;
Memo1.Lines.Add('Rel Error='+FloatToStr(Error));
finally
A.Free;
B.Free;
C.Free;
IBasis.Free;
Bi.Free;
end;
end;
|
|
Симплекс-метод, Минимизация целевой функции... 
Дата 24.3.2004, 18:29 (
Загрузка ...
