数学模型:    目标函数:Min S=C1X1+C2X2+....+CnXn
    约束条件:
    a11 X1 +a12 X2+...+a1n Xn≥(=,≤)b1
    a21 X1 +a22 X2+...+a2n Xn≥(=,≤)b2
    am1 X1 +am2 X2+...+amn Xn≥(=,≤)bm
    Xn >=0
Xn为决策变量,即各原料在配方中的比例;aij为技术系数,即各种原料相应的营养成分;bi为约束值,即配方应满足的各项营养指标或重量指标;Cj为成本系数,即每种原料的价格系数;m为约束条件的个数,n为变量个数。
 目的:使得在Xn有解的情况下,得到s.问题我简单的描述了一下,我想请你帮我写一个模块函数,用单纯形法即可,请指明输入输出的参数类型(要用到多维数组) ,避免使用指针,用object pascal写,我用的开发工具是delphi,希望能直接在程序中用,输出结果有两个  ,一个是总成本s,一个是数组参数Xn, 另外需要有一个判断是否有解的函数。

解决方案 »

  1.   

    有一个算法函数包,它是求标准线性模型的,不过它的目标函数是极大值,可以的话讲它转化一下是最好的了。{ @I 标示输入参数       @O  标示输出参数
    说明
    //M1 存放<=约束条件的个数  @I
    //M2 存放>=约束条件的个数  @I
    //M3 存放= 约束条件的个数  @I
    //M  约束条件的总个数  M=M1+M2+M3 @I
    //MP,NP 存储数组的物理维数,MP>=M+2,NP>=N+1
    //A (M+2)*(N+1)个元素的二维数组,其中A 的前M+1行是输入,输出参数,输入时 按
    // [  0   a01  a02      a0n
          b1  -a11  -a12    -a1n
          b2  -a21  -a22    -a2n      bm   -am1  -am2   -amn      bl    al1   al2    aln   ]IPOSV M个元素的一维整型数组,@O ,说明A的输出结果 若IZROV(J)=J,
    }
      

  2.   

    续://////////////////////////////////////////////////////////////////////////////////////////
    /////////////////////////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////////////////////////
    matrx2 = array of array of Extended;Procedure SIMPLX(var A:matrx2; M, N, MP, NP, M1, M2, M3:integer;
                          var ICASE:integer;var IZROV, IPOSV:array of integer);
    Label 1,2,3,4;
    const
        EPS  =0.000001;
    var
        L1,L2,L3:array[0..100] of integer;
        K,I,NL1,NL2,IR,KP,M12,IP,IQ,KH:integer;   Q1,BMAX:real;
    begin
        If M <> M1 + M2 + M3 Then
        begin
            ShowMessage(' Bad input constraint counts');
            Exit;
        end;
        NL1:=N;
        For K:=1 To N do
        begin
            L1[K]:=K;
            IZROV[K]:=K;
        end; 
        NL2:=M;
        For I:=1 To M do
        begin
            If A[I + 1, 1] < 0  Then
            begin
                ShowMessage('Bad input tableau.');
                Exit;
            End;
            L2[I]:=I;
            IPOSV[I]:=N + I;
        end;
        For I:=1 To M2 do
            L3[I]:=1;
        IR:=0;
        If M2 + M3 = 0 Then GoTo 3;
        IR:=1;
        For K:=1 To N + 1 do
        begin
            Q1:=0;
            For I:=M1 + 1 To M do
                Q1:=Q1 + A[I + 1, K];
            A[M + 2, K]:=-Q1;
        end;
        repeat
            SIMP1(A, MP, NP, M + 1, L1, NL1, 0, KP, BMAX);
            If (BMAX <= EPS) And (A[M + 2, 1] < -EPS) Then
            begin
              ICASE:=-1;
              Exit;
            end
            Else If (BMAX <= EPS) And (A[M + 2, 1] <= EPS) Then
            begin
              M12:=M1 + M2 + 1;
              If M12 <= M Then
              begin
                  For IP:=M12 To M do
                  begin
                      If IPOSV[IP] = IP + N Then
                      begin
                          SIMP1(A, MP, NP, IP, L1, NL1, 1, KP, BMAX);
                          If BMAX > 0  Then GoTo 1;
                      end;
                  end;
              end;
              IR:=0;
              M12:=M12 - 1;
              If M1 + 1 > M12 Then Exit;
              For I:=M1 + 1 To M12 do
              begin
                  If L3[I - M1] = 1 Then
                  begin
                      For K:=1 To N + 1 do
                          A[I + 1, K]:=-A[I + 1, K];
                  end;
              end;
              goto 3;
            end;
            SIMP2(A, M, N, MP, NP, L2, NL2, IP, KP, Q1);
            If IP = 0 Then
            begin
              ICASE:=-1;
              Exit;
            end;
      

  3.   

    续:1:      SIMP3(A, MP, NP, M + 1, N, IP, KP);
            If IPOSV[IP] >= N + M1 + M2 + 1 Then
            begin
              For K:=1 To NL1 do
                  If L1[K] = KP Then goto 4;
    4:        NL1:=NL1 - 1;
              For IQ:=K To NL1 do
                  L1[IQ]:=L1[IQ + 1];
            end
            Else
            begin
              If IPOSV[IP] < N + M1 + 1 Then GoTo 2;
              KH:=IPOSV[IP] - M1 - N;
              If L3[KH] = 0 Then GoTo 2;
              L3[KH]:=0;
            end;
            A[M + 2, KP + 1]:=A[M + 2, KP + 1] + 1;
            For I:=1 To M + 2 do
              A[I, KP + 1]:=-A[I, KP + 1];
    2:      IQ:=IZROV[KP];
            IZROV[KP]:=IPOSV[IP];
            IPOSV[IP]:=IQ;
        until IR <> 0;
    3:  SIMP1(A, MP, NP, 0, L1, NL1, 0, KP, BMAX);
        If BMAX <= 0  Then
        begin
            ICASE:=0;
            Exit;
        end;
        SIMP2(A, M, N, MP, NP, L2, NL2, IP, KP, Q1);
        If IP = 0 Then
        begin
            ICASE:=1;
            Exit;
        end;
        SIMP3(A, MP, NP, M, N, IP, KP);
        GoTo 2;
    end;
    Procedure SIMP1(var A:matrx2; MP, NP, MM:integer; LL:array of integer;
                              NLL, IABF:integer;var KP:integer; var BMAX:real);
    var
        K:integer;   TEST:real;
    begin
        KP:=LL[1];
        BMAX:=A[MM + 1, KP + 1];
        For K:=2 To NLL do
        begin
            If IABF = 0 Then
                TEST:=A[MM + 1, LL[K] + 1] - BMAX
            Else
                TEST:=Abs(A[MM + 1, LL[K] + 1]) - Abs(BMAX);
            If TEST > 0  Then
            begin
                BMAX:=A[MM + 1, LL[K] + 1];
                KP:=LL[K];
            end;
        end;
    end;
    Procedure SIMP2(var A:matrx2; M, N, MP, NP:integer; L2:array of integer;
                                     NL2:integer;var IP, KP:integer; Q1:real);
    label 1,2;
    const
        EPS = 0.000001;
    var
        I,J,II,K:integer;  FLAG,Q,QP,Q0:real;
    begin
        IP:=0;
        FLAG:=0;
        For I:=1 To NL2 do
        begin
            If A[L2[I] + 1, KP + 1] < -EPS Then FLAG:=1;
            If FLAG = 1 Then goto 1;
        end;
    1:  If FLAG = 0 Then Exit;
        Q1:=-A[L2[I] + 1, 1] / A[L2[I] + 1, KP + 1];
        IP:=L2[I];
        For I:=I + 1 To NL2 do
        begin
            II:=L2[I];
            If A[II + 1, KP + 1] < -EPS Then
            begin
                Q:=-A[II + 1, 1] / A[II + 1, KP + 1];
                If Q < Q1 Then
                begin
                    IP:=II;
                    Q1:=Q;
                end
                Else If Q = Q1 Then
                begin
                    For K:=1 To N do
                    begin
                        QP:=-A[IP + 1, K + 1] / A[IP + 1, KP + 1];
                        Q0:=-A[II + 1, K + 1] / A[II + 1, KP + 1];
                        If Q0 <> QP Then goto 2;
                    end;
    2:              If Q0 < QP Then IP:=II;
                end;
            end;
        end;
    end;Procedure SIMP3(var A:matrx2; MP, NP, I1, K1:integer;var IP, KP:integer);
    var
        PIV:real;  II,KK:integer;
    begin
        PIV:=1  / A[IP + 1, KP + 1];
        For II:=1 To I1 + 1 do
        begin
            If II - 1 <> IP Then
            begin
                A[II, KP + 1]:=A[II, KP + 1] * PIV;
                For KK:=1 To K1 + 1 do
                begin
                    If KK - 1 <> KP Then
                        A[II, KK]:=A[II, KK] - A[IP + 1, KK] * A[II, KP + 1];
                end;
            end;
        end;
        For KK:=1 To K1 + 1 do
            If KK - 1 <> KP Then A[IP + 1, KK]:=-A[IP + 1, KK] * PIV;
        A[IP + 1, KP + 1]:=PIV;
    end; 
     所有这些为验证上述子过程而编的验证过程
    implementation
    //PROGRAM D9R10
    //Driver for routine SIMPLX
    uses
      unit2;
      {$R *.DFM}procedure TForm1.Button1Click(Sender: TObject);
    const
      s1='%8.2f';
      N = 4; M = 4; NP = 5; MP = 6; M1 = 2; M2 = 1; M3 = 1; NM1M2 = N+M1+M2;
    var
      F:TextFile;   
      A:matrx2;
      IZROV,IPOSV:array[0..4] of integer;  
      ANUM:array[0..5] of real;
      TXT:array[0..7] of string;  
      ALPHA:array[0..5] of string;
      I,J,JJ,JMAX,ICASE:integer;  
      STR1,STR2,STR3:string;
    begin
      SetLength(A,7,6);
      TXT[1]:='x1';   TXT[2]:='x2';   TXT[3]:='x3';   TXT[4]:='x4';
      TXT[5]:='y1';   TXT[6]:='y2';   TXT[7]:='y3';
      A[1, 1]:=0;   A[1, 2]:=1;  A[1, 3]:=1;   A[1, 4]:=3;  A[1, 5]:=-0.5;
      A[2, 1]:=740; A[2, 2]:=-1; A[2, 3]:=0;   A[2, 4]:=-2; A[2, 5]:=0;
      A[3, 1]:=0;   A[3, 2]:=0;  A[3, 3]:=-2;  A[3, 4]:=0;  A[3, 5]:=7;
      A[4, 1]:=0.5; A[4, 2]:=0;  A[4, 3]:=-1;  A[4, 4]:=1;  A[4, 5]:=-2;
      A[5, 1]:=9;   A[5, 2]:=-1; A[5, 3]:=-1;  A[5, 4]:=-1; A[5, 5]:=-1;
      A[6, 1]:=0;   A[6, 2]:=0;  A[6, 3]:=0;   A[6, 4]:=0;  A[6, 5]:=0;
      SIMPLX(A, M, N, MP, NP, M1, M2, M3, ICASE, IZROV, IPOSV);
      //输出计算结果到文件
      AssignFile(F, 'd:\delphi_shu\p9\d9r10.dat');
      Rewrite(F);
      Writeln(F);
      If ICASE = 1 Then
        Writeln(F, '   Unbounded objective function')
      Else If ICASE = -1 Then
        Writeln(F, '   No solutions satisfy constraints given');
      JJ:=1;
      For I:=1 To N do
      begin
        If IZROV[I] <= N + M1 + M2 Then
        begin
          ALPHA[JJ]:=TXT[IZROV[I]];
          JJ:=JJ + 1;
        end;
      end;
      JMAX:=JJ - 1;
      Writeln(F);
      STR1:=' ';
      For JJ:=1 To JMAX do
        STR1:=STR1 + ALPHA[JJ]+'        ';
      Writeln(F, '                 ',STR1);
      For I:=1 To M + 1 do
      begin
        If I > 1 Then
          ALPHA[1]:=TXT[IPOSV[I - 1]]
        Else
          ALPHA[1]:='  ';
        ANUM[1]:=A[I, 1];
        JJ:=2;
        For J:=2 To N + 1 do
        begin
          If IZROV[J - 1] <= (N + M1 + M2) Then
          begin
            ANUM[JJ]:=A[I, J];
            JJ:=JJ + 1;
          end;
        end;
        JMAX:=JJ - 1;
        STR2:='  ';
        For JJ:=1 To JMAX do
        begin
          STR3:='   ';
          STR3:=COPY(FloatToStr(ANUM[JJ]),1,7);
          STR2:= STR2 + STR3 + '     ';
        end;
        Writeln(F, ALPHA[1],STR2);
      end;
      CloseFile(F);
      //屏幕显示计算结果
      memo1.Lines.LoadFromFile('d:\delphi_shu\p9\d9r10.dat');
    end;