Разработка программного модуля для нахождения оптимальных предельно-допустимых выбросов в атмосферу от группы источников

if (Simplex.Basis[i]<>-1) then Basis[i]:=Simplex.Basis[i]

else begin

SetAllLengths(N+1);

for j:=0 to M-1 do Cons[j].A[N-1]:=0;

Cons[i].A[N-1]:=1;

Cons[i].isT := true;

C[N-1] := 0;

for j:=0 to Simplex.N-1 do C[j] := C[j] + Simplex.Cons[i].A[j];

L := L + Cons[i].B;

end;

end;

end;

destructor TSimplex.Free;

begin

SetLength(C,0);

SetL

ength(Basis,0);

SetLength(Cons,0);

M:=0;

N:=0;

RealN := 0;

end;

function TSimplex.GetMin: extended;

var

i : integer;

begin

if (Max) then

Result := -L

else

Result := L;

end;

function TSimplex.GetSolution: TExtArray;

var

Solution : TExtArray;

i,j : integer;

begin

SetLength(Solution,RealN);

for j:=0 to RealN-1 do begin

Solution[j]:=0;

i:=0;

while ((i<M) and (Basis[i]<>j)) do inc(i);

if ((Basis[i]=j) and (i<M)) then Solution[j]:=Cons[i].B;

end;

Result:=Solution;

end;

procedure TSimplex.MulString(Number: integer; Value: extended);

var

j : integer;

begin

for j:=0 to N-1 do Cons[Number].A[j]:=Cons[Number].A[j]*Value;

Cons[Number].B:=Cons[Number].B*Value;

end;

procedure TSimplex.NormaliCe;

var

i : integer;

begin

for i:=0 to M-1 do if (Cons[i].Sign<>Equal) then begin

SetAllLengths(N+1);

if (Cons[i].Sign=Greater) then Cons[i].A[N-1]:=-1

else Cons[i].A[N-1]:=1;

Cons[i].Sign := Equal;

end;

end;

procedure TSimplex.SetAllLengths(Len: integer);

var

i, j : integer;

OldN : integer;

begin

OldN:=N;

N:=Len;

SetLength(C,N);

for i:=0 to M-1 do SetLength(Cons[i].A,N);

if (OldN<N) then begin

for j:=OldN to N-1 do begin

C[j]:=0;

for i:=0 to M-1 do Cons[i].A[j]:=0;

end;

end;

end;

function TSimplex.FoundInBasis(num:integer): integer;

var

i:integer;

f:boolean;

begin

f := false;

i := 0 ;

while (not f and (i<M)) do

begin

f := (Basis[i] = num);

inc(i);

end;

if (f) then

Result := i-1

else

Result := -1;

end;

function TSimplex.SimplexStep: integer;

var

i,j : integer;

f,opt : boolean;

x,y : integer; //координаты опорного элемента

CurMax : extended;

temp : array of TConstrain;

tempC : TExtArray;

begin

opt := true;

CurMax := -1;

for i := 0 to N-1 do

begin

//проверка на разрешимость

if (C[i] > 0) then

begin

opt := false; //а это попутная проверка на оптимальность

if (C[i] > CurMax) then //а это поиск ведущего столбца (максимальный элемент в C[i])

begin

CurMax := C[i];

x := i;

end;

f := true;

for j := 0 to M-1 do

f := f and (Cons[j].A[i] < 0);

if (f) then

begin

Result := SIMPLEX_NO_BOTTOM;

exit;

end;

end;

end;

if (opt) then

Result := SIMPLEX_DONE

else

begin

//зная номер ведущего столбца, ищем номер ведущей строки

CurMax := MaxExtended; //на самом деле тут будем искать минимум, а не Max

for j := 0 to M-1 do

if (Cons[j].A[x] > 0) then //идем только по положительным элементам

if (Cons[j].B/Cons[j].A[x] < CurMax) then

begin

CurMax := Cons[j].B/Cons[j].A[x];

y := j;

end

else if (DoPrec(Cons[j].B/Cons[j].A[x] - CurMax) = 0) then

if (Cons[j].isT) then

y := j;

//сохраняем текущие значения

SetLength(temp, M);

for j := 0 to M-1 do

begin

SetLength(temp[j].A, N);

for i := 0 to N-1 do

temp[j].A[i] := Cons[j].A[i];

temp[j].B := Cons[j].B;

end;

SetLength(tempC, N);

for i := 0 to N-1 do

tempC[i] := C[i];

//делаем пересчет таблицы

//строка делиться на ведущий элемент

MulString(y, 1/Cons[y].A[x]);

//преобразование остальных элементов

for j := 0 to M-1 do

begin

if (j <> y) then

begin

for i := 0 to N-1 do

begin

Cons[j].A[i] := DoPrec(temp[j].A[i] - temp[j].A[x]*temp[y].A[i]/temp[y].A[x]);

end;

Cons[j].B := DoPrec(temp[j].B - temp[j].A[x]*temp[y].B/temp[y].A[x]);

end

else

begin

for i := 0 to N-1 do

Cons[j].A[i] := DoPrec(Cons[j].A[i]);

end;

end;

//и строка с коэффициентами функции

for i := 0 to N-1 do

begin

C[i] := DoPrec(tempC[i] - tempC[x]*temp[y].A[i]/temp[y].A[x]);

end;

Basis[y] := x;

//и сама функция:

L := DoPrec(L - tempC[x]*temp[y].B/temp[y].A[x]);

for i:= 0 to M-1 do

SetLength(temp[i].A, 0);

SetLength(temp, 0);

SetLength(tempC, 0);

Result := SIMPLEX_NEXT_STEP;

end;

end;

function TSimplex.Solve: integer;

var

i,j : integer;

Simplex : TSimplex;

f : boolean;

Step : integer;

cc : extended;

begin

//oldN := N;

NormaliCe;

f:=false;

if (not CheckBasis) then begin

Simplex:=TSimplex.CreateBasis(self);

Simplex.Solve;

f:=Simplex.GetMin<>0;

if (not f) then for i:=0 to M-1 do begin

for j:=0 to N-1 do Cons[i].A[j]:=Simplex.Cons[i].A[j];

Cons[i].B:=Simplex.Cons[i].B;

Cons[i].isT := false;

Basis[i]:=Simplex.Basis[i];

cc := C[Basis[i]];

for j:=0 to N-1 do

C[j] := DoPrec(C[j] - cc*Cons[i].A[j]);

L := DoPrec(L - cc*Cons[i].B);

end;

Simplex.Free;

end;

if (f) then Step:=SIMPLEX_NO_SOLUTION

else repeat

Step:=SimplexStep;

until (Step<>SIMPLEX_NEXT_STEP);

//SetAllLengths(OldN);

Result:=Step;

end;

{ TIntSimplex }

constructor TIntSimplex.Create(_C:TExtArray; MaximiCe:boolean=false);

begin

//CurFound:=false;

inherited;

end;

function TIntSimplex.GetIntMin: extended;

begin

Result:=GetMin;

end;

function TIntSimplex.GetIntSolution: TExtArray;

begin

Result:=GetSolution;

end;

function TIntSimplex.IsInteger(Value:extended):boolean;

begin

Result:=((Value=floor(Value)) or (Value=ceil(Value)));

end;

function TIntSimplex.IntSolve: integer;

var

i : integer;

OldN : integer;

FractCol : integer;

FractRow : integer;

TmpX : TExtArray;

TmpCons : TExtArray;

NewValue : extended;

begin

if (Solve=SIMPLEX_DONE) then begin

//if (not CurFound or ((Simplex.GetMin<CurL) and not Max) or ((Simplex.GetMin>CurL) and Max)) then begin

TmpX:=GetSolution;

i:=0;

while ((i<RealN) and IsInteger(TmpX[i])) do inc(i);

FractCol:=i;

if (FractCol<>RealN) then begin // если найдена хотя бы одна нецелая переменная

OldN:=N;

SetLength(TmpCons,N);

FractRow := FoundInBasis(FractCol);

for i := 0 to N-1 do

if (FoundInBasis(i) = -1) then

TmpCons[i] := Cons[FractRow].A[i] - Floor(Cons[FractRow].A[i])

else

TmpCons[i] := 0;

NewValue := Cons[FractRow].B - Floor(Cons[FractRow].B);

//if (Max) then

AddCons(NewValue, TmpCons, Greater);

//else

// AddCons(NewValue, TmpCons, Less);

Result := IntSolve;

SetAllLengths(OldN); // удаляем пустые столбцы в конце, если они есть

 Скачать реферат