unit1.cpp

lp求解的单纯形算法

           i+=1;
           if(text[i]=='=')  i+=1;
           while(text[i]!='-'&&text[i]!='+'&&!IsNumber(text[i]))
                 {
                    if(IsLetter(text[i]))
                       {
                         MessageDlg("对不起!输入错误!约束不等式右端不能有变量!",
                                    mbOK,TMsgDlgButtons()<<mbOK,0);
                         return;
                       }
                    i+=1;
                 }
           sign='+';
           if(text[i]=='-')
                  sign='-';
           while(!IsNumber(text[i])&&i<total)  i+=1;
           s=0;
           while((IsNumber(text[i])||text[i]=='.')&&i<total)
                   {  NumberWork[s]=text[i]; s+=1; i+=1; }
           NumberWork[s]='\0';
           RHS[p]= StrToFloat(NumberWork);
           if(sign=='-')   RHS[p]*=-1;
        }
    //约束条件已经处理完毕
    delete [] text;

//******************************************************************************
   //最大转化为最小
    if(type)
        for(i=0;i<Col;i++)
                 C[i]*=-1;

    //把Name复制给AnsiName,方便后面的操作
    Ansi_Vector AnsiName;
    AnsiName.Length=Col+2*Row;
    for(i=0;i<Col;i++)
        AnsiName[i]=Name[i];

    //AnsiName已经代替Name记录变量名,Name需要释放
    for(i=0;i<Col;i++)
           free(Name[i]);
    free(Name);

    //扩展成标准问题
    k=0;
    for(i=0;i<Row;i++)
       {
         if(signs[i]=='>')
            {  A[i][Col+k]=-1;
               AnsiName[Col+k]="Std"+String(k+1);
               k+=1;
            }
         else if(signs[i]=='<')
            {  A[i][Col+k]=1;
               AnsiName[Col+k]="Std"+String(k+1);
               k+=1;
            }
       }

    int ColStd=Col+k;//标准问题的约束矩阵列数
    int ColAssist=ColStd+Row;//辅助问题的约束矩阵的列数
    
    //扩展成辅助问题,为第一阶段准备
    for(i=0;i<Row;i++)
       {
          A[i].Length-=(Row-k);
          A[i][ColStd+i]=1;
          AnsiName[ColStd+i]="Assist"+String(i+1);//人工变量
          if(RHS[i]<0)//确保右端向量RHS总是大于零
              {
                 RHS[i]*=-1;
                 for(j=0;j<ColStd;j++)
                      A[i][j]*=-1;
              }
       }
    C.Length-=(Row-k);

//******************************************************************************
//下面开始实行两阶段法的第一阶段来寻找标准问题的基本可行基，以便第二阶段的进行
    Int_Vector BaseRecord;//用来记录基变量在AnsiName中的下标
    BaseRecord.Length=Row;//基矩阵是Row*Row阵
    //辅助问题初始基变量是Row个人工变量,它们在Name最后Row个位置
    for(i=0;i<Row;i++)
           BaseRecord[i]=ColStd+i;

    Double_Vector AssistTest,StdTest;
    AssistTest.Length=ColAssist+1;
    StdTest.Length=ColAssist+1;//这里最后一个分量是用来存放最优值的
    //以下求出辅助问题的初始检验向量
    for(i=0;i<ColStd;i++)
        {
           AssistTest[i]=0;
           for(j=0;j<Row;j++)
               AssistTest[i]+=A[j][i];
           StdTest[i]=-1*C[i];
        }
    for(i=ColStd;i<=ColAssist;i++)  {  AssistTest[i]=0; StdTest[i]=0;  }
    for(i=0;i<Row;i++)    AssistTest[ColAssist]+=RHS[i];

    double max,min,temp;
    int BaseIn,BaseOut,ResultType,Step=0;//入基,出基,结果类型

    while(1)
         {
             Step+=1;
             StatusBar1->Panels->Items[6]->Text="正在寻找基本可行解.....";
             BaseIn=Max_Index(AssistTest,ColAssist);//AssistTest前ColAssist分量
                                                   //最大值的下标
             max=AssistTest[BaseIn];//入基

             if(max<=Difference)//已经找到最优解
                   break;

             max=Max(A,BaseIn,Row);//A的k列前Row个分量的最大值

             if(max<=Difference)
                   break;//无界,辅助问题中无此情况

             min=Infinite;
             for(i=0;i<Row;i++)
                 if(A[i][BaseIn]>Difference)
                     if(RHS[i]/A[i][BaseIn]< min)//出基
                         {  BaseOut=i;  min= RHS[i]/ A[i][BaseIn]; }

             temp=A[BaseOut][BaseIn]; //以 A[BaseOut][BaseIn]为转轴元
             for(j=0;j<ColAssist;j++) //进行旋转变换
                     A[BaseOut][j]/=temp;
             RHS[BaseOut]/=temp;

             for(i=0;i<Row;i++)
                  if(i!=BaseOut)
                     {
                        temp=A[i][BaseIn];
                        for(j=0;j<ColAssist;j++)
                             A[i][j]-=(A[BaseOut][j]*temp);
                        RHS[i]-=(RHS[BaseOut]*temp);
                     }

             temp=AssistTest[BaseIn];
             for(j=0;j<ColAssist;j++)
                  AssistTest[j]-=(A[BaseOut][j]*temp);
             AssistTest[ColAssist]-=(RHS[BaseOut]*temp);

             temp=StdTest[BaseIn];
             for(j=0;j<ColAssist;j++)
                  StdTest[j]-=(A[BaseOut][j]*temp);
             StdTest[ColAssist]-=(RHS[BaseOut]*temp);

             BaseRecord[BaseOut]=BaseIn;//换基
         }
         
    double g=AssistTest[ColAssist];//辅助问题的最优值
    if(fabs(g)<=Difference)    AssistTest[ColAssist]=0;//消除计算误差
    else
        {
           StatusBar1->Panels->Items[6]->Text="没有基本可行解....已停止....";
           MessageDlg("原问题没有基本可行解!", mbOK,TMsgDlgButtons()<<mbOK,0);
           return;
        }
    //辅助问题的最优值是零
    k=-1;
    for(i=0;i<Row;i++)//判断基变量中是否有人工变量
        if(BaseRecord[i]>=ColStd)
            { k=i;  break;  }

    while(k!=-1)//基变量中有人工变量
        {

           StatusBar1->Panels->Items[6]->Text="正在消去人工变量....";
           s=-1;
           for(j=0;j<ColStd;j++)
              if(  fabs(A[k][j])>Difference )
                       {  s=j; break;  }

           if(s!=-1)//k行元素不全为零,则要以A[k][s]为转轴元进行旋转变换
               {
                  Step+=1;
                  BaseOut=k;  BaseIn=s;//此处必有RHS[BaseOut]==0
                  temp=A[BaseOut][BaseIn];
                  for(j=0;j<ColAssist;j++)
                        A[BaseOut][j]/=temp;

                  for(i=0;i<Row;i++)
                       if(i!=BaseOut)
                          {
                             temp=A[i][BaseIn];
                             for(j=0;j<ColAssist;j++)
                                  A[i][j]-=(A[BaseOut][j]*temp);
                          }

                  temp=AssistTest[BaseIn];
                  for(j=0;j<ColAssist;j++)
                       AssistTest[j]-=(A[BaseOut][j]*temp);

                  temp=StdTest[BaseIn];
                  for(j=0;j<ColAssist;j++)
                       StdTest[j]-=(A[BaseOut][j]*temp);

                  BaseRecord[BaseOut]=BaseIn;
               }

           else//k行元素全为零,有多余的约束行,应删去
               {
                  for(i=k;i<Row-1;i++)
                      {
                         RHS[i]=RHS[i+1];
                         BaseRecord[i]=BaseRecord[i+1];
                         for(j=0;j<ColAssist;j++)
                               A[i][j]=A[i+1][j];
                      }
                  Row=Row-1;
               }

           k=-1;
           for(i=0;i<Row;i++)//判断基变量中是否有人工变量
                  if(BaseRecord[i]>=ColStd)
                         { k=i;  break;  }
        }
    //基变量中已经没有人工变量了
//******************************************************************************
//下面开始实行两阶段法的第二阶段
    while(1)
         {
              Step+=1;
              StatusBar1->Panels->Items[6]->Text="正在寻找最优解.....";
              BaseIn=Max_Index(StdTest,ColStd);
              max=StdTest[BaseIn];

              if(max<=Difference)//已经找到最优解
                { ResultType=1;  break; }

              max=Max(A,BaseIn,Row);//A的k列前Row个分量的最大值

              if(max<=Difference)
                { ResultType=0;  break; }//无界

              min=Infinite;
              for(i=0;i<Row;i++)
                 if(A[i][BaseIn]>Difference)
                     if(RHS[i]/A[i][BaseIn]< min)//出基
                         {  BaseOut=i;  min= RHS[i]/ A[i][BaseIn]; }

              temp=A[BaseOut][BaseIn]; //以 A[BaseOut][BaseIn]为转轴元
              for(j=0;j<ColStd;j++) //进行旋转变换
                     A[BaseOut][j]/=temp;
              RHS[BaseOut]/=temp;

              for(i=0;i<Row;i++)
                  if(i!=BaseOut)
                     {
                        temp=A[i][BaseIn];
                        for(j=0;j<ColStd;j++)
                             A[i][j]-=(A[BaseOut][j]*temp);
                        RHS[i]-=(RHS[BaseOut]*temp);
                     }

              temp=StdTest[BaseIn];
              for(j=0;j<ColStd;j++)
                     StdTest[j]-=(A[BaseOut][j]*temp);
              StdTest[ColAssist]-=(RHS[BaseOut]*temp);

              BaseRecord[BaseOut]=BaseIn;
         }
    if(ResultType==0)
       {
            StatusBar1->Panels->Items[6]->Text="原问题无界.....已停止....";
            MessageDlg("原问题无界!", mbOK,TMsgDlgButtons()<<mbOK,0);
            return;
       }
    else
       {
         StatusBar1->Panels->Items[6]->Text="已经找到最优解,停止." ;
         Double_Vector x;
         double z=StdTest[ColAssist];
         if(type)  z*=-1;
         x.Length=Col;
         for(i=0;i<Col;i++)  x[i]=0;
         for(i=0;i<Row;i++)  x[BaseRecord[i]]=RHS[i];
         TResultForm *Result;
         if(NowChild->Result==NULL)
               NowChild->Result=new TResultForm(Application);
         Result=NowChild->Result;
         Result->MemoRe->Lines->Clear();
         Result->MemoRe->Lines->Add(" 你输入问题经过旋转次数为:");
         Result->MemoRe->Lines->Add("                         "+IntToStr(Step));
         Result->MemoRe->Lines->Add(" 最优值是： ");
         Result->MemoRe->Lines->Add("         "+ FloatToStr(z) );
         for(i=0;i<Col;i++)
             Result->MemoRe->Lines->Add(" "+AnsiName[i]+
                                           "**************"+FloatToStr(x[i]) );
       }
}


//--------------------------------------------帮助文档--------------------------
void __fastcall TMainForm::HelpClick(TObject *Sender)
{
    if(HaveHelp==false)
      {
        THelpForm *help;
        help=new THelpForm(Application);
        HaveHelp=true;
      }
}

//-----------------程序退出-----------------------------------------------------


void __fastcall TMainForm::ExitClick(TObject *Sender)
{
     Close();
}
//--------------------------------OVER------------------------------------------