单纯形法性规划问题求解器 -- 算法驿站-- 编程爱好者博客

我花了一天时间写的，为了做作业，不用用笔在纸上算了。简单写了一个分析字符串的东西，没有专门的考虑词法分析什么的，只从中间提取变量名和参数表，有一个符号表功能，矩阵运算是原来数值分析实习时写的CMatrix。 LPSolver 单纯形法线性规划问题求解器 v0.1 rickone [实现功能]1,直接单纯形法2,大M法,M近似3,两阶段法 [输入格式]MAX或MIN(目标函数=目标函数表达式){ 约束方程1; 约束方程2; ...}1,目标函数变量为第一个扫描到的变量标识 2,线性表达式的每一项:<系数><变量名> 3,约束形式:<= = >==不能写成== 4,符号'='为方程左右的边界标识符 5,符号';'为方程结束标识符,最后一个也不能省 6,隐含的约束条件有,所有的变量>=0 [例子] MAX(z=x1+x2+x3){ x1-x2<=100; 2x1-1.5x2>=10; x1+3x3=50;} QQ:85104865email to: rickone@sina.comdate:2007/03/19 VC6源码（部分，仅含主窗口的代码，分析字符串和LP算法，Matrix的只给个接口说明吧）： //: Matrix.h#if !defined(AFX_MATRIX_H__83023DF6_28C2_4C65_A12D_9DAAF6624E84__INCLUDED_)#define AFX_MATRIX_H__83023DF6_28C2_4C65_A12D_9DAAF6624E84__INCLUDED_ #if _MSC_VER > 1000#pragma once#endif // _MSC_VER > 1000// Matrix.h : header file// /////////////////////////////////////////////////////////////////////////////// CMatrix window class CMatrix : public CObject{// Support Serialization DECLARE_SERIAL(CMatrix)// Constructionpublic: CMatrix(); CMatrix(int n); CMatrix(int height,int width); CMatrix(CMatrix&); CMatrix(double* da,int height,int width); CMatrix(CMatrix& mx,int ys,int xs,int height,int width); // Attributespublic: double *data; int w,h; // Operations CMatrix& operator=(const CMatrix &m);public: // Overrides // ClassWizard generated virtual function overrides //{{AFX_VIRTUAL(CMatrix) //}}AFX_VIRTUAL void Serialize(CArchive& ar); // Implementationpublic: virtual ~CMatrix(); enum{LIN=0,COL=1};//行和列的标识 void size(int& height,int& width); double& elem(int y,int x);//选取元素[i,j]，引用返回 void display(); CString GetMatrix(); void swap(int t1,int t2,int tag=LIN);//交换两行或两列，默认为行变换 void mul(int t1,double c,int tag=LIN);//某一行(列)乘上c倍 void pt(int t1,int t2,double c,int tag=LIN);//初等行或列变化，默认为行变换 void unin(const CMatrix&,int tag=LIN); // Generated message map functionsprotected: //{{AFX_MSG(CMatrix) // NOTE - the ClassWizard will add and remove member functions here. //}}AFX_MSG}; ///////////////////////////////////////////////////////////////////////////// //{{AFX_INSERT_LOCATION}}// Microsoft Visual C++ will insert additional declarations immediately before the previous line. #endif // !defined(AFX_MATRIX_H__83023DF6_28C2_4C65_A12D_9DAAF6624E84__INCLUDED_) //: LPSolverDlg.cpp// LPSolverDlg.cpp : implementation file// #include "stdafx.h"#include "LPSolver.h"#include "LPSolverDlg.h"#include "DataInputDlg.h"#include "ArgTable.h"#include "math.h" #ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif /////////////////////////////////////////////////////////////////////////////// CAboutDlg dialog used for App About class CAboutDlg : public CDialog{public: CAboutDlg(); // Dialog Data //{{AFX_DATA(CAboutDlg) enum { IDD = IDD_ABOUTBOX }; //}}AFX_DATA // ClassWizard generated virtual function overrides //{{AFX_VIRTUAL(CAboutDlg) protected: virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV support //}}AFX_VIRTUAL // Implementationprotected: //{{AFX_MSG(CAboutDlg) //}}AFX_MSG DECLARE_MESSAGE_MAP()}; CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD){ //{{AFX_DATA_INIT(CAboutDlg) //}}AFX_DATA_INIT} void CAboutDlg::DoDataExchange(CDataExchange* pDX){ CDialog::DoDataExchange(pDX); //{{AFX_DATA_MAP(CAboutDlg) //}}AFX_DATA_MAP} BEGIN_MESSAGE_MAP(CAboutDlg, CDialog) //{{AFX_MSG_MAP(CAboutDlg) // No message handlers //}}AFX_MSG_MAPEND_MESSAGE_MAP() /////////////////////////////////////////////////////////////////////////////// CLPSolverDlg dialog CLPSolverDlg::CLPSolverDlg(CWnd* pParent /*=NULL*/) : CDialog(CLPSolverDlg::IDD, pParent){ //{{AFX_DATA_INIT(CLPSolverDlg) m_strShow = _T("LPSolver 0.1 rickone\r\n"); m_bShowStep = TRUE; m_strBigM = _T("1000"); m_nMethod = 1; //}}AFX_DATA_INIT // Note that LoadIcon does not require a subsequent DestroyIcon in Win32 m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME); m_strProblem = _T("");} void CLPSolverDlg::DoDataExchange(CDataExchange* pDX){ CDialog::DoDataExchange(pDX); //{{AFX_DATA_MAP(CLPSolverDlg) DDX_Control(pDX, IDC_EDIT_SHOW, m_edit); DDX_Control(pDX, IDC_LIST_ARG, m_listbox); DDX_Text(pDX, IDC_EDIT_SHOW, m_strShow); DDX_Check(pDX, IDC_SHOW_STEP, m_bShowStep); DDX_Text(pDX, IDC_EDIT_BIGM, m_strBigM); DDX_Radio(pDX, IDC_RADIO1, m_nMethod); //}}AFX_DATA_MAP} BEGIN_MESSAGE_MAP(CLPSolverDlg, CDialog) //{{AFX_MSG_MAP(CLPSolverDlg) ON_WM_SYSCOMMAND() ON_WM_PAINT() ON_WM_QUERYDRAGICON() ON_BN_CLICKED(IDC_DATA_INPUT, OnDataInput) ON_BN_CLICKED(IDC_CLOSE_WINDOW, OnCloseWindow) ON_BN_CLICKED(IDC_CLEAN_TEXT, OnCleanText) ON_WM_CTLCOLOR() //}}AFX_MSG_MAPEND_MESSAGE_MAP() /////////////////////////////////////////////////////////////////////////////// CLPSolverDlg message handlers BOOL CLPSolverDlg::OnInitDialog(){ CDialog::OnInitDialog(); // Add "About..." menu item to system menu. // IDM_ABOUTBOX must be in the system command range. ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX); ASSERT(IDM_ABOUTBOX < 0xF000); CMenu* pSysMenu = GetSystemMenu(FALSE); if (pSysMenu != NULL) { CString strAboutMenu; strAboutMenu.LoadString(IDS_ABOUTBOX); if (!strAboutMenu.IsEmpty()) { pSysMenu->AppendMenu(MF_SEPARATOR); pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu); } } // Set the icon for this dialog. The framework does this automatically // when the application's main window is not a dialog SetIcon(m_hIcon, TRUE); // Set big icon SetIcon(m_hIcon, FALSE); // Set small icon // TODO: Add extra initialization here return TRUE; // return TRUE unless you set the focus to a control} void CLPSolverDlg::OnSysCommand(UINT nID, LPARAM lParam){ if ((nID & 0xFFF0) == IDM_ABOUTBOX) { CAboutDlg dlgAbout; dlgAbout.DoModal(); } else { CDialog::OnSysCommand(nID, lParam); }} // If you add a minimize button to your dialog, you will need the code below// to draw the icon. For MFC applications using the document/view model,// this is automatically done for you by the framework. void CLPSolverDlg::OnPaint() { if (IsIconic()) { CPaintDC dc(this); // device context for painting SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0); // Center icon in client rectangle int cxIcon = GetSystemMetrics(SM_CXICON); int cyIcon = GetSystemMetrics(SM_CYICON); CRect rect; GetClientRect(&rect); int x = (rect.Width() - cxIcon + 1) / 2; int y = (rect.Height() - cyIcon + 1) / 2; // Draw the icon dc.DrawIcon(x, y, m_hIcon); } else { CDialog::OnPaint(); }} // The system calls this to obtain the cursor to display while the user drags// the minimized window.HCURSOR CLPSolverDlg::OnQueryDragIcon(){ return (HCURSOR) m_hIcon;} void CLPSolverDlg::OnOK() { // TODO: Add extra validation here UpdateData(); GetMatrixFromStr(m_strProblem); ShowArgTable(); if(m_matrix.h>0 && m_matrix.w>0) { BOOL bRun=FALSE; if(m_nMethod==0) { double BM=atof(m_strBigM); if(BM<50) m_strShow+=">!大M值太小了<<\r\n"; else bRun=BigMSolver(m_matrix,BM,m_bShowStep); } else if(m_nMethod==1) { bRun=TwoStepSolver(m_matrix,m_bShowStep); } else if(m_nMethod==2) { bRun=NormalSolver(m_matrix,m_bShowStep); } if(bRun) { CString tmp; tmp.Format(">目标函数最优解为 %lf\r\n",((CArgElem*)m_table.m_arr.GetAt(0))->value); m_strShow+=">计算结果见变量列表<<\r\n"+tmp; } ShowArgTable(); } else { m_strShow+=">!单纯形表未准备好<<\r\n"; } UpdateData(FALSE);} void CLPSolverDlg::OnDataInput() { // TODO: Add your control notification handler code here CDataInputDlg dlg; if(!(m_strProblem=="")) dlg.m_strSource=m_strProblem; if(IDOK==dlg.DoModal()) { m_strProblem=dlg.m_strSource; GetMatrixFromStr(dlg.m_strSource); ShowArgTable(); m_strShow+="--LP Problem Math Description--\r\n"+dlg.m_strSource+"\r\n--end--\r\n"+m_matrix.GetMatrix(); UpdateData(FALSE); }} void CLPSolverDlg::OnCleanText() { // TODO: Add your control notification handler code here UpdateData(); m_strShow=_T("LPSolver 0.1 rickone\r\n"); UpdateData(FALSE);} void CLPSolverDlg::OnCloseWindow() { // TODO: Add your control notification handler code here EndDialog(IDOK);} void CLPSolverDlg::Analyze(CString &strSrc){ /* int nSize=m_arrayArgTable.GetSize(); for(int i=0;i<nSize;++i) { delete m_arrayArgTable.GetAt(i); } m_arrayArgTable.RemoveAll(); int index=0;//变元编号 //寻找MIN or MAX int pos=strSrc.Find("MAX"); int max=1; if(pos==-1) { pos=strSrc.Find("MIN"); max=-1; } if(pos==-1) { return; } //目标函数 pos=strSrc.Find('(',pos); int endpos=strSrc.Find('=',pos); CString ss=strSrc.Mid(pos+1,endpos-pos-1); m_arrayArgTable.Add(new CArgTable(ss,0.0,index++));*/} //DEL void CLPSolverDlg::DivArgVal(CString str, CString &val, CString &nam)//DEL {//DEL int nLen=str.GetLength();//DEL int i;//DEL for(i=0;i<nLen;++i)//DEL {//DEL char c=str.GetAt(i);//DEL if((c<'0' || c>'9') && c!='.')//DEL break;//DEL }//DEL val=str.Left(i);//DEL nam=str.Mid(i);//DEL } void CLPSolverDlg::ShowArgTable(){ m_listbox.ResetContent(); for(int i=0;i<m_table.n;++i) { CArgElem *p=(CArgElem*)(m_table.m_arr.GetAt(i)); CString tmp; tmp.Format("%-10s%-10.3lf%d",p->name,p->value,p->index); m_listbox.InsertString(m_listbox.GetCount(),tmp); }} int CLPSolverDlg::GetArgFromStr(CString &strSrc){ m_table.CleanTable(); int pos; int rindex=0; int fnum=0; pos=strSrc.Find("MIN"); m_bMin=TRUE; if(pos==-1) { pos=strSrc.Find("MAX"); m_bMin=FALSE; if(pos==-1) return -1; } CString tmp; int nLen=strSrc.GetLength(); for(int i=pos+3;i<nLen;++i) { int j,cc,nSize; CString arg; char c=strSrc.GetAt(i); switch(c) { case '(': case ')': case '{': case '}': case '+': case '-': case '<': case '=': case '>': case ';': case '\r': case '\n': tmp.TrimLeft(); tmp.TrimRight(); nSize=tmp.GetLength(); for(j=0;j<nSize;++j) { cc=tmp.GetAt(j); if((cc<'0' || cc>'9') && cc!='.') break; } arg=tmp.Mid(j); if(arg.IsEmpty()) { tmp.Empty(); break; } if(!m_table.IsExist(arg)) m_table.InsertArg(arg); tmp.Empty(); break; default: tmp+=c; } if(c=='}') break; else if(c=='<' || c=='>') { //加入松驰变量 arg.Format("RX%d",rindex++); m_table.InsertArg(arg); } else if(c==';') { //标识方程个数 fnum++; } } return fnum;} void CLPSolverDlg::GetMatrixFromStr(CString &strSrc){ int h=GetArgFromStr(strSrc)+1; int w=m_table.n; CMatrix matrix(h,w); // 不含人工变量的单纯形表 , 初始全0 int pos; int rindex=0; double mark=1.0,side=1.0; pos=strSrc.Find("MIN"); if(pos==-1) { pos=strSrc.Find("MAX"); side=-1.0; if(pos==-1) return; } CString tmp; CString val,arg; int nLen=strSrc.GetLength(); int line=0; for(int i=pos+3;i<nLen;++i) { int j,cc,nSize; CString arg; char c=strSrc.GetAt(i); switch(c) { case '(': case ')': case '{': case '}': case '+': case '-': case '<': case '=': case '>': case ';': case '\r': case '\n': tmp.TrimLeft(); tmp.TrimRight(); nSize=tmp.GetLength(); for(j=0;j<nSize;++j) { cc=tmp.GetAt(j); if((cc<'0' || cc>'9') && cc!='.') break; } val=tmp.Left(j); arg=tmp.Mid(j); if(val.IsEmpty() && arg.IsEmpty()) { tmp.Empty(); break; } double valf; if(val.IsEmpty()) valf=1.0; else valf=atof(val); int aindex; if(arg.IsEmpty()) aindex=0; // 常数列 else { aindex=m_table.GetArgIndex(arg); if(aindex==0) { //目标函数z valf=0.0; } } matrix.elem(line,aindex)=mark*side*valf; mark=1.0; tmp.Empty(); break; default: tmp+=c; } switch(c) { case '}': i=nLen; break; case '<': case '>': //加入松驰变量 arg.Format("RX%d",rindex++); matrix.elem(line,m_table.GetArgIndex(arg))=(c=='<'?1.0:-1.0); break; case ';': case ')': case '=': if(c!='=') { //标识方程个数 line++; side=1.0; } else { //转换到=另一边 side=-side; } break; case '-': //负号 mark=-1.0; break; } } m_matrix=matrix;} BOOL CLPSolverDlg::BigMSolver(CMatrix &matrix, double M, BOOL bShowStep){ //const double NegBigM=-500; //增加人工变量 CString arg; int aIndex=0; CMatrix hm(matrix.h,matrix.h-1); m_strShow+=""; for(int i=0;i<hm.w;++i) { hm.elem(0,i)=-M; hm.elem(i+1,i)=1.0; arg.Format("HU%d",aIndex++); m_table.InsertArg(arg); } CMatrix mx=matrix; mx.unin(hm); if(bShowStep) { arg.Format(">大M法计算模块\r\nM=%lf 近似\r\n--LPSolver--\r\n",M); m_strShow+=arg+">增加人工变量后>>\r\n"+mx.GetMatrix(); } return NormalSolver(mx,bShowStep);} BOOL CLPSolverDlg::TwoStepSolver(CMatrix &matrix, BOOL bShowStep){ //第一阶段 double *z=new double[matrix.w]; int i; for(i=0;i<matrix.w;++i) { z[i]=matrix.elem(0,i); // 保存原目标函数 matrix.elem(0,i)=0.0; } CString arg; int aIndex=0; CMatrix hm(matrix.h,matrix.h-1); for(i=0;i<hm.w;++i) { hm.elem(0,i)=-1.0; hm.elem(i+1,i)=1.0; arg.Format("HU%d",aIndex++); m_table.InsertArg(arg); } CMatrix mx1=matrix; mx1.unin(hm); if(bShowStep) { m_strShow+=">两阶段法计算模块\r\n>增加人工变量后>>\r\n" +mx1.GetMatrix()+">第一阶段计算开始>>\r\n"; } NormalSolver(mx1,bShowStep); // 解第一阶段 if(fabs(mx1.elem(0,0))<1e-30) { //第二阶段 CMatrix mx2(mx1,0,0,matrix.h,matrix.w); for(i=0;i<mx2.w;++i) { mx2.elem(0,i)=z[i]; } if(bShowStep) { m_strShow+=">第二阶段计算开始>>\r\n>去掉人工变量后>>\r\n"+mx2.GetMatrix(); } return NormalSolver(mx2,bShowStep); } // 解不为0,无可行解 if(bShowStep) { m_strShow+=">!第一阶段的解不为0，原LP问题无解<<\r\n"; } delete[] z; return FALSE;} BOOL CLPSolverDlg::NormalSolver(CMatrix &matrix, BOOL bShowStep){ // 直接对单纯形表进行求解 int *BaseLin=new int[matrix.h]; BOOL *BaseCol=new BOOL[matrix.w]; double *phi=new double[matrix.h]; CString tmp; for(int j=0;j<matrix.w;++j) { BaseCol[j]=FALSE; } for(int i=1;i<matrix.h;++i) { for(int j=1;j<matrix.w;++j) { if(fabs(matrix.elem(i,j)-1.0)<1e-30) { for(int k=1;k<matrix.h && (k==i?1:fabs(matrix.elem(k,j))<1e-30);++k); if(k==matrix.h) { //找到一个单位基向量 matrix.pt(i,0,-matrix.elem(0,j)); // 消去基向量从目标函数 BaseCol[j]=TRUE; BaseLin[i]=j; break; } } } if(j==matrix.w) { if(bShowStep) { m_strShow+=">!找不到初始基可行解,算法无法启动<<\r\n"; } //无初始基可行解 delete[]BaseCol; delete[]BaseLin; delete[]phi; return FALSE; } } if(bShowStep) { m_strShow+=">从目标函数中消去基变量>>\r\n"+matrix.GetMatrix(); } while(TRUE) { double max=-1; int maxp; for(int i=1;i<matrix.w;++i) { if(!BaseCol[i] && matrix.elem(0,i)>=max) { max=matrix.elem(0,i); maxp=i; } } if(max<=0.0)//非基变量参数全为负，找到最优解 break; if(bShowStep) { tmp.Format(">找到换入变量%s(第%d列)\r\n",m_table.GetArgName(maxp),maxp+1); m_strShow+=tmp; } //否则计算phi值，同时用最小原则找minphi值 double minph=1e100; int minphp=-1; for(int j=1;j<matrix.h;++j) { if(matrix.elem(j,maxp)>0.0) { phi[j]=-matrix.elem(j,0)/matrix.elem(j,maxp); if(phi[j]<minph) { minph=phi[j]; minphp=j; } } } if(minphp==-1) { if(bShowStep) { m_strShow+=">!解无界<<\r\n"; } //无界解 delete[]BaseCol; delete[]BaseLin; delete[]phi; return FALSE; } if(bShowStep) { tmp.Format(">找到换出变量%s(第%d列)\r\n",m_table.GetArgName(BaseLin[minphp]),BaseLin[minphp]+1); m_strShow+=tmp; } BaseCol[BaseLin[minphp]]=FALSE; //换出 BaseLin[minphp]=maxp; //换入 //用初等行变换将(minphp,maxp)消元成1 matrix.mul(minphp,1.0/matrix.elem(minphp,maxp)); for(int k=0;k<matrix.h;++k) { if(k!=minphp) { matrix.pt(minphp,k,-matrix.elem(k,maxp)); } } if(bShowStep) { m_strShow+=">初等行变换后的结果>>\r\n"+matrix.GetMatrix(); } } if(bShowStep) { m_strShow+=">完成<<\r\n"; } // 将结果显示在变量列表中 double fResult; if(m_bMin) fResult=-matrix.elem(0,0); else fResult=matrix.elem(0,0); ((CArgElem*)m_table.m_arr.GetAt(0))->value=fResult; for(i=1;i<matrix.h;++i) { ((CArgElem*)m_table.m_arr.GetAt(BaseLin[i]))->value=-matrix.elem(i,0); } delete[]BaseCol; delete[]BaseLin; delete[]phi; return TRUE;} HBRUSH CLPSolverDlg::OnCtlColor(CDC* pDC, CWnd* pWnd, UINT nCtlColor) { HBRUSH hbr = CDialog::OnCtlColor(pDC, pWnd, nCtlColor); // TODO: Change any attributes of the DC here if(pWnd->GetDlgCtrlID()==IDC_EDIT_SHOW) { hbr=(HBRUSH)GetStockObject(BLACK_BRUSH); pDC->SetBkColor(RGB(0,0,0)); pDC->SetTextColor(RGB(0,128,0)); } // TODO: Return a different brush if the default is not desired return hbr;}

博客介绍

正文

单纯形法性规划问题求解器2007-03-19 22:09:00

评论