数模论坛

 找回密码
 注-册-帐-号
搜索
热搜: 活动 交友 discuz
查看: 4148|回复: 8

本人前几天封装的单纯性法解线性规划的类

[复制链接]
发表于 2004-4-12 22:22:12 | 显示全部楼层 |阅读模式
/*
  Name:                                LP.h
  Copyright: Knight-Studio
  Author: Ran Ni
  Date: 12-04-04 14:30
  Description: The Declaration of the Class LinearProgramming
*/
#include <vector>
#ifndef _LINEARPROGRAMMING_H_
#define _LINEARPROGRAMMING_H_

using namespace std;

class LinearProgramming
{
  public:
    LinearProgramming(vector<vector<double> > a,vector<double> bb,vector<double> cc):A(a),b(bb),c(cc),Feasible(false),BestSolution(bb){}
    vector<double> GetBestSolution()
    {
           return BestSolution;
    }
    double GetMini()
    {
           return Minimum;
    }
    bool Run();   
        
  protected:
    vector<vector<double> >A;
    vector<double> b;
    vector<double> c;
    double Minimum;
    bool Feasible;
    vector<double> BestSolution;
};
#endif


 楼主| 发表于 2004-4-12 22:23:05 | 显示全部楼层
/*
  Name:                                  LP.cpp
  Copyright: Knight-Studio
  Author: Ran Ni
  Date: 12-04-04 14:29
  Description: The Definition of Class LinearProgramming
*/
#include <vector>
#include <math.h>
#include <algorithm>
#include "LP.h"

using namespace std;

void Unit(vector<double> &x)
{
        for(int i=0;i<(int)x.size();++i)
        {
               x=0;
        }
}
      
double operator * (vector<double> a,vector<double> b)
{
        double data(0);
        for(int i=0;i<(int)a.size();++i)
        {
               data+=a*b;
        }
        return data;
}

vector<vector<double> > Inverse(vector<vector<double> > a)
{
        vector<vector<double> > b(a.size());
        for(int i=0;i<(int)b.size();++i)
        {
               vector<double> c(a.size(),0);
               c=1;
               b=c;
        }
        for(int k=0;k<(int)b.size();++k)
        {
               int g=k;
               for(int i=k;i<(int)b.size();++i)
               {
                      if(fabs(a[k])>fabs(a[k][g]))
                      {
                             g=i;
                      }
               }
               if(g!=k)
               {
                      for(int i=0;i<(int)b.size();++i)
                      {
                             double temp=a[k];
                             a[k]=a[g];
                             a[g]=temp;
                             temp=b[k];
                             b[k]=b[g];
                             b[g]=temp;
                      }
               }
               if(a[k][k]!=0)
               {
                      double x=a[k][k];
                      for(int i=0;i<(int)b.size();++i)
                      {
                             a[k]=(a[k]/x);
                             b[k]=(b[k]/x);
                      }
                      for(int i=0;i<(int)b.size();++i)
                      {
                             if(i!=k)
                             {
                                   double y=a[k];
                                   for(int j=k;j<(int)b.size();++j)
                                   {
                                         a[j]=a[j]-a[j][k]*y/a[k][k];                                         
                                   }
                                   for(int j=0;j<(int)b.size();++j)
                                   {
                                         b[j]=b[j]-b[j][k]*y/a[k][k];                                       
                                   }
                             }
                      }
               }
        }
        return b;
}
               
vector<double> operator * (vector<vector<double> > a,vector<double> b)
{
        vector<double> x(b.size(),0);
        for(int i=0;i<(int)x.size();++i)
        {
               double data(0);
               for(int j=0;j<(int)b.size();++j)
               {
                      data+=b[j]*a[j];
               }
               x=data;
        }
        return x;
}

vector<double> CalculateW(vector<double> c,vector<vector<double> > b)
{
       vector<double> data(b.size(),0);
       for(int i=0;i<(int)data.size();++i)
       {
               data=c*b;
       }
       return data;
}

vector<double> operator - (vector<double> a,vector<double> b)
{
       for(int i=0;i<(int)a.size();++i)
       {
               a-=b;
       }
       return a;
}

int max(vector<double> c)
{
       int k(0);
       for(int i=0;i<(int)c.size();++i)
       {
               if(c>c[k])
               {
                       k=i;
               }
       }
       return k;
}

bool Judge(vector<double> a)
{
       int i=0;
       while((a<=0)&&(i<(int)a.size()))
       {
               i++;
       }
       if(i==(int)a.size())
       {
               return true;
       }
       else
       {
               return false;
       }
}

bool LinearProgramming::Run()
{
       vector<vector<double> > AA(A);
       for(int i=0;i<(int)b.size();++i)
       {
               vector<double> tem(b.size(),0);
               tem=1;
               AA.push_back(tem);
       }
       vector<double> cc(c);
       for(int i=0;i<(int)b.size();++i)
       {
               cc.push_back(0);
       }
       vector<int> Base(b.size(),0);
       for(int i=0;i<(int)b.size();++i)
       {
               Base=A.size()+i;
       }
       vector<double> cb(b.size(),0);
       vector<vector<double> > Br;
       vector<double> xb;
       vector<double> xx(AA.size(),0);
       double f;
       vector<double> w;
       vector<double> dd(AA.size(),0);
step1: vector<vector<double> > B;
       for(int i=0;i<(int)b.size();++i)
       {
               B.push_back(AA[Base]);
       }      
       for(int i=0;i<(int)b.size();++i)
       {
               cb=(cc[Base]);
       }
       Br=Inverse(B);
       xb=(Br*b);
       Unit(xx);
       for(int i=0;i<(int)Base.size();++i)
       {
               xx[Base]=xb;
       }
       f=cb*xb;
       w=(CalculateW(cb,Br));
       Unit(dd);
       for(int i=0;i<AA.size();++i)
       {
              if(find(Base.begin(),Base.end(),i)==Base.end())
              {
                      dd=(w*AA)-cc;
              }
       }
       int k(max(dd));
       if(dd[k]<=0)
       {
              this->Feasible=true;
              BestSolution=xx;
              Minimum=f;
              goto step2;
       }
       else
       {
              vector<double> yk(Br*AA[k]);
              if(Judge(yk))
              {
                     this->Feasible=false;
                     goto step2;              
              }
              else
              {
                     int r(0);
                     while(!(yk[r]>0))
                     {
                            ++r;
                     }
                     double dddd(b[r]/yk[r]);
                     for(int i=r;i<(int)yk.size();++i)
                     {
                            if((yk>0)&&((b/yk)<dddd))
                            {
                                    r=i;
                                    dddd=b/yk;
                            }
                     }
                     Base[r]=k;
                     goto step1;
              }
       }
step2: if(this->Feasible==true)
       {
              return true;
       }
       else
       {
              return false;
       }
}

 楼主| 发表于 2004-4-12 22:24:12 | 显示全部楼层
程序不是很长,因为用了很多的容器变量和范型算法。
发表于 2004-4-13 02:14:17 | 显示全部楼层
顶!!!
发表于 2004-4-13 19:51:08 | 显示全部楼层
赞一个
发表于 2004-4-25 05:44:58 | 显示全部楼层
<>好厉害  我是新手  多指教!</P>
发表于 2004-4-25 06:39:00 | 显示全部楼层
<>帮knight代发他改进后的代码:</P><>/*
  Name: LinearProgramming
  Copyright: Knight-Studio
  Author: Ran Ni
  Date: 12-04-04 14:29
  Description: The Definition of Class LinearProgramming
*/
#include &lt;vector&gt;
#include &lt;math.h&gt;
#include &lt;algorithm&gt;
#include "LP.h"</P><>using namespace std;</P><P>void Unit(vector&lt;double&gt; &amp;x)
{
        for(int i=0;i&lt;(int)x.size();++i)
        {
               x<i>=0;
        }
}
      
double operator * (vector&lt;double&gt; a,vector&lt;double&gt; b)
{
        double data(0);
        for(int i=0;i&lt;(int)a.size();++i)
        {
               data+=a<i>*b<i>;
        }
        return data;
}</P><P>vector&lt;vector&lt;double&gt; &gt; Inverse(vector&lt;vector&lt;double&gt; &gt; a)
{
        vector&lt;vector&lt;double&gt; &gt; b(a.size());
        for(int i=0;i&lt;(int)b.size();++i)
        {
               vector&lt;double&gt; c(a.size(),0);
               c<i>=1;
               b<i>=c;
        }
        for(int k=0;k&lt;(int)b.size();++k)
        {
               int g=k;
               for(int i=k;i&lt;(int)b.size();++i)
               {
                      if(fabs(a[k]<i>)&gt;fabs(a[k][g]))
                      {
                             g=i;
                      }
               }
               if(g!=k)
               {
                      for(int i=0;i&lt;(int)b.size();++i)
                      {
                             double temp=a<i>[k];
                             a<i>[k]=a<i>[g];
                             a<i>[g]=temp;
                             temp=b<i>[k];
                             b<i>[k]=b<i>[g];
                             b<i>[g]=temp;
                      }
               }
               if(a[k][k]!=0)
               {
                      double x=a[k][k];
                      for(int i=0;i&lt;(int)b.size();++i)
                      {
                             a<i>[k]=(a<i>[k]/x);
                             b<i>[k]=(b<i>[k]/x);
                      }
                      for(int i=0;i&lt;(int)b.size();++i)
                      {
                             if(i!=k)
                             {
                                   double y=a[k]<i>;
                                   for(int j=k;j&lt;(int)b.size();++j)
                                   {
                                         a[j]<i>=a[j]<i>-a[j][k]*y/a[k][k];                                         
                                   }
                                   for(int j=0;j&lt;(int)b.size();++j)
                                   {
                                         b[j]<i>=b[j]<i>-b[j][k]*y/a[k][k];                                       
                                   }
                             }
                      }
               }
        }
        return b;
}
               
vector&lt;double&gt; operator * (vector&lt;vector&lt;double&gt; &gt; a,vector&lt;double&gt; b)
{
        vector&lt;double&gt; x(b.size(),0);
        for(int i=0;i&lt;(int)x.size();++i)
        {
               double data(0);
               for(int j=0;j&lt;(int)b.size();++j)
               {
                      data+=b[j]*a[j]<i>;
               }
               x<i>=data;
        }
        return x;
}</P><P>vector&lt;double&gt; CalculateW(vector&lt;double&gt; c,vector&lt;vector&lt;double&gt; &gt; b)
{
       vector&lt;double&gt; data(b.size(),0);
       for(int i=0;i&lt;(int)data.size();++i)
       {
               data<i>=c*b<i>;
       }
       return data;
}</P><P>vector&lt;double&gt; operator - (vector&lt;double&gt; a,vector&lt;double&gt; b)
{
       for(int i=0;i&lt;(int)a.size();++i)
       {
               a<i>-=b<i>;
       }
       return a;
}</P><P>int max(vector&lt;double&gt; c)
{
       int k(0);
       for(int i=0;i&lt;(int)c.size();++i)
       {
               if(c<i>&gt;c[k])
               {
                       k=i;
               }
       }
       return k;
}</P><P>int min(vector&lt;double&gt; c)
{
       int k(0);
       for(int i=0;i&lt;(int)c.size();++i)
       {
               if(c<i>&lt;c[k])
               {
                       k=i;
               }
       }
       return k;
}</P><P>bool Judge(vector&lt;double&gt; a)
{
       int i=0;
       while((a<i>&lt;=0)&amp;&amp;(i&lt;(int)a.size()))
       {
               i++;
       }
       if(i==(int)a.size())
       {
               return true;
       }
       else
       {
               return false;
       }
}</P><P>bool Jugdeb(vector&lt;double&gt; b)
{
       int i=0;
       while((b<i>&gt;=0)&amp;&amp;(i&lt;(int)b.size()))
       {
               i++;
       }
       if(i==(int)b.size())
       {
               return true;
       }
       else
       {
               return false;
       }
}</P><P>bool LinearProgramming::Run()
{
       vector&lt;vector&lt;double&gt; &gt; AA(A);
       Num=0;
       for(int i=0;i&lt;(int)b.size();++i)
       {
               vector&lt;double&gt; tem(b.size(),0);
               tem<i>=1;
               AA.push_back(tem);
       }
       vector&lt;double&gt; cc(c);
       int k;
       for(int i=0;i&lt;(int)b.size();++i)
       {
               cc.push_back(0);
       }
       vector&lt;int&gt; Base(b.size(),0);
       for(int i=0;i&lt;(int)b.size();++i)
       {
               Base<i>=A.size()+i;
       }
       vector&lt;double&gt; cb(b.size(),0);
       vector&lt;vector&lt;double&gt; &gt; Br;
       vector&lt;double&gt; xb;
       vector&lt;double&gt; xx(AA.size(),0);
       double f;
       vector&lt;double&gt; w;
       vector&lt;double&gt; xt(AA.size()+1,0);
       vector&lt;double&gt; dd(AA.size(),0);
       vector&lt;double&gt; bb(b);
       vector&lt;vector&lt;double&gt; &gt; B(b.size());
       if(!Jugdeb(bb))
       {
                vector&lt;double&gt; xa(b.size(),-1);
                AA.push_back(xa);
                vector&lt;double&gt; ct(cc.size(),0);
                ct.push_back(1);
                int rt(min(bb));
                for(int i=0;i&lt;(int)AA.size();++i)
                {
                         AA<i>[rt]=-AA<i>[rt];
                }
                bb[rt]=-bb[rt];
                for(int i=0;i&lt;(int)bb.size();++i)
                {
                         if(i!=rt)
                         {
                                  for(int j=0;j&lt;(int)AA.size();++j)
                                  {
                                           AA[j]<i>-=(AA[j][rt]*(AA[AA.size()-1]<i>/AA[AA.size()-1][rt]));
                                  }
                                  bb<i>+=(bb[rt]);
                         }
                }
                Base[rt]=AA.size()-1;
                ++Num;
step3:          vector&lt;vector&lt;double&gt; &gt; Bs;
                for(int i=0;i&lt;(int)Base.size();++i)
                {
                         Bs.push_back(AA[Base<i>]);
                }      
                for(int i=0;i&lt;(int)bb.size();++i)
                {
                         cb<i>=(ct[Base<i>]);
                }
                Br=Inverse(Bs);
                xb=(Br*bb);
                Unit(xt);
                for(int i=0;i&lt;(int)Base.size();++i)
                {
                         xt[Base<i>]=xb<i>;
                }
                w=(CalculateW(cb,Br));
                Unit(dd);
                for(int i=0;i&lt;(int)AA.size();++i)
                {
                         if(find(Base.begin(),Base.end(),i)==Base.end())
                         {
                                  dd<i>=(w*AA<i>)-ct<i>;
                         }
                }
                int ks(max(dd));
                if(dd[ks]&lt;=0)
                {
                        goto step1;
                }
                else
                {
                        vector&lt;double&gt; yk(Br*AA[ks]);
                        int r(0);
                        while(!(yk[r]&gt;0))
                        {
                                ++r;
                        }
                        double dddd(bb[r]/yk[r]);
                        for(int i=r;i&lt;(int)yk.size();++i)
                        {
                                if((yk<i>&gt;0)&amp;&amp;((bb<i>/yk<i>)&lt;dddd))
                                {
                                         r=i;
                                         dddd=bb<i>/yk<i>;
                                }
                        }
                        Base[r]=ks;
                        ++Num;
                        goto step3;                        
                }               
       }
step1: if(find(Base.begin(),Base.end(),AA.size()-1)!=Base.end())
       {
               if(xt[AA.size()-1]!=0)
               {
                        Feasible=false;
                        goto step2;
               }
               else
               {
                        vector&lt;int&gt; :: iterator it=find(Base.begin(),Base.end(),(int)AA.size()-1);
                        int i=0;
                        while(Base<i>!=(*it))
                        {
                                ++i;
                        }
                        int j=0;
                        while(!((AA[j]<i>!=0)&amp;&amp;(find(Base.begin(),Base.end(),j)==Base.end())))
                        {
                                ++j;
                        }
                        *it=j;
               }
       }
       for(int i=0;i&lt;(int)B.size();++i)
       {
               B<i>=(AA[Base<i>]);
       }      
       for(int i=0;i&lt;(int)bb.size();++i)
       {
               cb<i>=(cc[Base<i>]);
       }
       Br=Inverse(B);
       xb=(Br*bb);
       Unit(xx);
       for(int i=0;i&lt;(int)Base.size();++i)
       {
               xx[Base<i>]=xb<i>;
       }
       f=cb*xb;
       w=(CalculateW(cb,Br));
       Unit(dd);
       for(int i=0;i&lt;(int)AA.size()-1;++i)
       {
              if(find(Base.begin(),Base.end(),i)==Base.end())
              {
                     dd<i>=(w*AA<i>)-cc<i>;
              }
       }
       k=max(dd);
       if(dd[k]&lt;=0)
       {
              this-&gt;Feasible=true;
              this-&gt;Limited=true;
              BestSolution=xx;
              Minimum=f;
              goto step2;
       }
       else
       {
              vector&lt;double&gt; yk(Br*AA[k]);
              if(Judge(yk))
              {
                     this-&gt;Feasible=false;
                     this-&gt;Limited=false;
                     goto step2;              
              }
              else
              {
                     int r(0);
                     while(!(yk[r]&gt;0))
                     {
                            ++r;
                     }
                     double dddd(bb[r]/yk[r]);
                     for(int i=r;i&lt;(int)yk.size();++i)
                     {
                            if((yk<i>&gt;0)&amp;&amp;((bb<i>/yk<i>)&lt;dddd))
                            {
                                    r=i;
                                    dddd=bb<i>/yk<i>;
                            }
                     }
                     Base[r]=k;
                     ++Num;
                     goto step1;
              }
       }
step2: if(this-&gt;Feasible==true)
       {
              return true;
       }
       else
       {
              return false;
       }
}
</P>
发表于 2004-4-25 06:39:32 | 显示全部楼层
<>/*
  Name: LinearPrograming
  Copyright: Knight-Studio
  Author: Ran Ni
  Date: 12-04-04 14:30
  Description: The Declaration of the Class LinearProgramming
*/
#include &lt;vector&gt;
#ifndef _LINEARPROGRAMMING_H_
#define _LINEARPROGRAMMING_H_</P><>using namespace std;</P><>class LinearProgramming
{
  public:
    LinearProgramming(vector&lt;vector&lt;double&gt; &gt; a,vector&lt;double&gt; bb,vector&lt;double&gt; cc):A(a),b(bb),c(cc),Feasible(false),BestSolution(bb),Limited(true),Num(0){}
    vector&lt;double&gt; GetBestSolution()
    {
           return BestSolution;
    }
    double GetMini()
    {
           return Minimum;
    }
    bool GetLimited()
    {
           return Limited;
    }
    int GetNum()
    {
           return Num;
    }
    bool Run();   
        
  protected:
    vector&lt;vector&lt;double&gt; &gt;A;
    vector&lt;double&gt; b;
    vector&lt;double&gt; c;
    bool Limited;
    double Minimum;
    bool Feasible;
    int Num;
    vector&lt;double&gt; BestSolution;
};
#endif</P>
发表于 2005-8-22 00:49:53 | 显示全部楼层
<>强烈支持!</P>
您需要登录后才可以回帖 登录 | 注-册-帐-号

本版积分规则

小黑屋|手机版|Archiver|数学建模网 ( 湘ICP备11011602号 )

GMT+8, 2024-11-27 10:44 , Processed in 0.059288 second(s), 19 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表