前言
由于列车运行图问题涉及整数分支定界算法,而该算法需要单纯形法计算每个分支最优解,故此处运用C#实现单纯形法。希望大佬批评指点。
说明
Winform窗体,只有一个Button按钮;其中,第一部分代码为算法实现部分,第二部分代码为调用部分,调用部分c为原问题目标函数系数,a为约束系数,b为资源系数,d较为特殊,本文规定,若为<=约束,d值为-1,若为=约束,d值为0,若为>=约束,d值为1。适用问题为最小值问题,若为最大值,请将目标函数乘以-1即可。可直接粘贴调试。
第一部分
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace SimplexMethod
{
/// <summary>
/// 单纯形算法
/// </summary>
class SimplexMethodAlgorithm
{
double[] CMat;//系数矩阵
double[][] AMat;//约束矩阵
double[] bMat;//资源矩阵
int MConstraintCont;//需要添加人工变量的约束个数
double M = 5000000;//大M法中的M值
int[] XBIndex;//定义存储基变量标号的一维数组
int[] XNIndex;//定义存储非基变量标号的一维数组
double[] XBValue;//定义基变量值的一维数组
double[] XNEnter;//定义进基变量的判别数
double[] XBExit;//定义出基变量的判别数
double[] CBValue;//定义基变量的目标函数系数的向量
double[] CNValue;//定义非基变量的目标函数系数的向量
double solutionF = double.MaxValue;//目标函数值
//double[] CBValueMultiplicationInverseofB;//定义CB与B-1的乘积值
//double[] InverseofBMultiplicationP;//定义B-1与P的乘积值
//double[,] B;//定义存储基变量在约束表达式中的系数矩阵
int t;//定义当前循环次数
int enterIndex;//定义当前进基序号
int exitIndex;//定义当前出基序号
//int[] exitIndexJudgeSolutionDegeneracy;//定义当前出基序号及判断是否退化
int[] judgeSymbolDirection;//判断约束的符号方向,"<="为-1,“=”为0,“>="为1
public List<double> solutionSet = new List<double>();//存储最优解列表集
bool isContinue = true;
/// <summary>
/// 算法执行函数
/// </summary>
public void ExecuteAlgorithm()
{
Console.WriteLine("\r\n\r\n\r\n");
Console.WriteLine("算法开始迭代:");
for (t = 1; t < 10; t++)//第t次迭代
{
Console.WriteLine("第 " + t + " 次迭代!");
GetXNEnter(CBValue, XNIndex);//获取非基变量检验数
GetEnterIndex(XNEnter);//计算进基变量下标
if (enterIndex != -1)//当算法需继续迭代时
{
GetXBExit(enterIndex);//获取基变量检验数
if (isContinue)//重新计算AMat、bMat
{
GetExitIndex(XBExit);//计算出基变量下标
//判断是否存在退化情况
List<int> exitMinIndex = new List<int>(); //存储最小比值下标
for (int i = 0; i < XBExit.Length; i++)
{
if (XBExit[i] == XBExit.Min())
{
exitMinIndex.Add(XBIndex[i]);
}
}
if (exitMinIndex.Count > 1)//如果出现退化情况--运用Bland法则重新选择出入基变量元素
{
for (int j = 0; j < XBIndex.Length; j++)//查找索引下标最小元素为出基变量
{
if (XBIndex[j] == exitMinIndex.Min())
{
exitIndex = j;
break;
}
}
for (int j = 0; j < XNEnter.Length; j++)//查找索引下标最小元素为入基变量
{
if (XNEnter[j] < 0)
{
enterIndex = XNIndex[j];
break;
}
}
}
UpdateAmatandbMatData();//更新A、b数据
XBIndex[exitIndex] = enterIndex;//更新基变量下标数组
CBValue[exitIndex] = CMat[enterIndex];//更新基变量对应系数值
int xnIndex = 0;
for (int i = 0; i < CMat.Length; i++)//更新非基变量下标数组
{
if (!XBIndex.Contains(i))
{
XNIndex[xnIndex++] = i;
}
}
for (int i = 0; i < XNIndex.Length; i++)//更新CNValue数组
{
CNValue[i] = CMat[XNIndex[i]];
}
for (int i = 0; i < bMat.Length; i++)//更新基变量值数组
{
XBValue[i] = bMat[i];
}
}
else
{
break;
}
}
else
{
break;
}
solutionSet.Add(solutionF);
solutionF = 0;
for (int i = 0; i < CBValue.Length; i++)//计算迭代解
{
solutionF += CBValue[i] * bMat[i];
}
}
}
/// <summary>
/// 更新约束矩阵AMat、资源矩阵bMat
/// </summary>
private void UpdateAmatandbMatData()
{
double douValue = AMat[exitIndex][enterIndex];//定位exitIndex行XNIndex[enterIndex]列的值
for (int i = 0; i < CMat.Length; i++)//第exitIndex行AMat值线性变换
{
AMat[exitIndex][i] /= douValue;
}
double[][] AMatCopy = new double[bMat.Length][];//此处须存储约束矩阵AMat
for (int i = 0; i < AMatCopy.GetLength(0); i++)
{
AMatCopy[i] = new double[CMat.Length];
for (int j = 0; j < AMatCopy[i].Length; j++)
{
AMatCopy[i][j] = AMat[i][j];
}
}
bMat[exitIndex] /= douValue;//第exitIndex行bMat值线性变换
for (int i = 0; i < bMat.Length; i++)//除exitIndex行之外的AMat值线性变换
{
if (i != exitIndex)
{
for (int j = 0; j < CMat.Length; j++)//更新每个约束系数值
{
AMat[i][j] += AMatCopy[exitIndex][j] * (-AMatCopy[i][enterIndex]);
}
bMat[i] += bMat[exitIndex] * (-AMatCopy[i][enterIndex]);
}
}
}
/// <summary>
/// 计算出基变量索引
/// </summary>
/// <param name="_xbExit">出基变量检验数</param>
private void GetExitIndex(double[] _xbExit)
{
int minIndex = 0;
double douMinValue = _xbExit[0];
for (int i = 1; i < bMat.Length; i++)
{
if (douMinValue >= _xbExit[i])
{
douMinValue = _xbExit[i];
minIndex = i;
}
}
exitIndex = minIndex;
}
/// <summary>
/// 获取出基变量检验数
/// </summary>
/// <param name="_enterIndex">进基变量索引值</param>
private void GetXBExit(int _enterIndex)
{
int judgeExistSolution = 0;
for (int i = 0; i < bMat.Length; i++)//判断进基约束列向量中小于等于0的元素个数
{
if (AMat[i][_enterIndex] <= 0)
{
judgeExistSolution++;
}
}
if (judgeExistSolution == bMat.Length)//如果所有元素均小于等于0
{
Console.WriteLine("**********该问题无可行解!**********");
exitIndex = -1;
isContinue = false;
return;
}
for (int i = 0; i < bMat.Length; i++)
{
if (AMat[i][_enterIndex] <= 0)
{
XBExit[i] = double.MaxValue;
}
else
{
XBExit[i] = bMat[i] / AMat[i][_enterIndex];
}
isContinue = true;
}
}
/// <summary>
/// 计算进基变量索引
/// </summary>
/// <param name="_xnEnter">非基变量检验数</param>
private void GetEnterIndex(double[] _xnEnter)
{
List<int> minute = new List<int>();//负值检验数标号集合----可进基标号集合
List<int> positive = new List<int>();//正值检验数标号集合
List<int> zero = new List<int>();//零值检验数标号集合
for (int i = 0; i < _xnEnter.Length; i++)
{
if (_xnEnter[i] < 0)
{
minute.Add(i);
}
else if (_xnEnter[i] == 0)
{
zero.Add(i);
}
else
{
positive.Add(i);
}
}
if (minute.Count == 0)//如果非基变量检验数均大于等于0时,则表示当前解就是最优解或无解
{
for (int j = 0; j < CBValue.Length; j++)
{
if (CBValue[j] == M)
{
if (XBValue[j] != 0)
{
enterIndex = -1;
}
}
}
if (positive.Count == _xnEnter.Length)//如果当前解为最优解
{
for (int i = 0; i < bMat.Length; i++)
{
Console.Write("X" + (XBIndex[i] + 1).ToString());
Console.Write(" ");
}
Console.WriteLine();
for (int i = 0; i < bMat.Length; i++)
{
Console.Write(bMat[i].ToString());
Console.Write(" ");
}
Console.WriteLine();
Console.Write(solutionF);
Console.WriteLine();
Console.WriteLine("**********找到最优解!**********");
enterIndex = -1;
}
else
{
for (int i = 0; i < bMat.Length; i++)
{
Console.Write("X" + (XBIndex[i] + 1).ToString());
Console.Write(" ");
}
Console.WriteLine();
for (int i = 0; i < bMat.Length; i++)
{
Console.Write(bMat[i].ToString());
Console.Write(" ");
}
Console.WriteLine();
Console.Write(solutionF);
Console.WriteLine();
Console.WriteLine("**********有无穷最优解!**********");
enterIndex = -1;
}
}
else//如果非基变量检验数存在负值,该解还不是最优解,需要继续迭代,选择检验数最小
{
double minValue = _xnEnter[0];
int minValueIndex = 0;
for (int j = 1; j < _xnEnter.Length; j++)
{
if (minValue > _xnEnter[j])
{
minValueIndex = j;
minValue = _xnEnter[j];
}
}
enterIndex = XNIndex[minValueIndex];
}
}
/// <summary>
/// 计算非基变量检验数
/// </summary>
/// <param name="_cbvalue">基变量目标函数系数向量</param>
/// <param name="_xnindex">非基变量标号向量</param>
private void GetXNEnter(double[] _cbvalue, int[] _xnindex)
{
for (int i = 0; i < _xnindex.Length; i++)//针对每个非基变量资源向量值
{
XNEnter[i] = 0;
for (int j = 0; j < _cbvalue.Length; j++)//针对每个基变量系数值
{
XNEnter[i] += _cbvalue[j] * AMat[j][_xnindex[i]];
}
XNEnter[i] = CMat[_xnindex[i]] - XNEnter[i];
}
}
/// <summary>
/// 构造函数
/// </summary>
///<param name="_cMat">目标函数系数向量</param>
///<param name="_aMat">约束方程系数矩阵</param>
///<param name="_bMat">约束方程资源向量</param>
///<param name="_judgeSymbolDirection">约束条件的符号方向表示向量</param>
public SimplexMethodAlgorithm(double[] _cMat, double[][] _aMat, double[] _bMat, int[] _judgeSymbolDirection)
{
CMat = _cMat;
AMat = _aMat;
bMat = _bMat;
judgeSymbolDirection = _judgeSymbolDirection;
XBIndex = new int[bMat.Length];
//OutParamater(CMat, AMat, bMat,XBIndex);
Console.WriteLine("\r\n\r\n\r\n");
GetStandardEquation(ref CMat, ref AMat, ref judgeSymbolDirection);//标准化方程
OutParamater(CMat, AMat, bMat, XBIndex);//输出标准化方程
XNIndex = new int[CMat.Length - XBIndex.Length];
int xnIndex = 0;
for (int i = 0; i < CMat.Length; i++)
{
if (!XBIndex.Contains(i))
{
XNIndex[xnIndex++] = i;
}
}
XBValue = new double[XBIndex.Length];
XNEnter = new double[XNIndex.Length];
XBExit = new double[XBIndex.Length];
CBValue = new double[XBIndex.Length];
CNValue = new double[XNIndex.Length];
solutionF = 0;
//CBValueMultiplicationInverseofB = new double[CBValue.Length];
//InverseofBMultiplicationP = new double[CBValue.Length];
for (int i = 0; i < XBValue.Length; i++)
{
XBValue[i] = bMat[i];
CBValue[i] = CMat[XBIndex[i]];
}
for (int i = 0; i < XNIndex.Length; i++)
{
CNValue[i] = CMat[XNIndex[i]];
}
for (int i = 0; i < CBValue.Length; i++)
{
solutionF += CBValue[i] * bMat[i];
}
solutionSet.Add(solutionF);
t = 0;
enterIndex = -1;
exitIndex = -1;
//exitIndexJudgeSolutionDegeneracy = new int[2];
}
/// <summary>
/// 方程标准化
/// </summary>
/// <param name="_cMat">目标函数系数向量</param>
/// <param name="_aMat">约束条件系数矩阵</param>
/// <param name="_judgeSymbolDirection">约束条件符号表示向量</param>
private void GetStandardEquation(ref double[] _cMat, ref double[][] _aMat, ref int[] _judgeSymbolDirection)
{
int index = 0;//基向量索引下标
for (int i = 0; i < _judgeSymbolDirection.Length; i++)
{
//小于等于0的约束
if (_judgeSymbolDirection[i] == -1)
{
//更新A[i]值
double[][] ARowCopy = new double[_aMat.GetLength(0)][];//定义拷贝并拓展A值多维数组
for (int j = 0; j < ARowCopy.GetLength(0); j++)
{
ARowCopy[j] = new double[_cMat.Length + 1];
if (i == j)
{
for (int k = 0; k < _cMat.Length; k++)
{
ARowCopy[j][k] = _aMat[j][k];
}
ARowCopy[j][_cMat.Length] = 1;
}
else
{
for (int k = 0; k < _cMat.Length; k++)
{
ARowCopy[j][k] = _aMat[j][k];
}
}
}
_aMat = ARowCopy;
//更新C值
double[] CRowCopy = new double[_cMat.Length + 1];
for (int j = 0; j < CRowCopy.Length; j++)
{
if (j < _cMat.Length)
{
CRowCopy[j] = _cMat[j];
}
else
{
CRowCopy[j] = 0;
}
}
_cMat = CRowCopy;
//存储基变量下标
XBIndex[index++] = _cMat.Length - 1;
}
//等于0的约束
else if (_judgeSymbolDirection[i] == 0)
{
//更新A[i]值
double[][] ARowCopy = new double[_aMat.GetLength(0)][];//定义拷贝并拓展A值多维数组
for (int j = 0; j < ARowCopy.GetLength(0); j++)
{
ARowCopy[j] = new double[_cMat.Length + 1];
if (i == j)
{
for (int k = 0; k < _cMat.Length; k++)
{
ARowCopy[j][k] = _aMat[j][k];
}
ARowCopy[j][_cMat.Length] = 1;
}
else
{
for (int k = 0; k < _cMat.Length; k++)
{
ARowCopy[j][k] = _aMat[j][k];
}
}
}
_aMat = ARowCopy;
//更新C值
double[] CRowCopy = new double[_cMat.Length + 1];
for (int j = 0; j < CRowCopy.Length; j++)
{
if (j < _cMat.Length)
{
CRowCopy[j] = _cMat[j];
}
else
{
CRowCopy[j] = M;
}
}
_cMat = CRowCopy;
//存储基变量下标
XBIndex[index++] = _cMat.Length - 1;
}
//大于等于0的约束
else
{
//更新A[i]值
double[][] ARowCopy = new double[_aMat.GetLength(0)][];//定义拷贝并拓展A值多维数组
for (int j = 0; j < ARowCopy.GetLength(0); j++)
{
ARowCopy[j] = new double[_cMat.Length + 2];
if (i == j)
{
for (int k = 0; k < _cMat.Length; k++)
{
ARowCopy[j][k] = _aMat[j][k];
}
ARowCopy[j][_cMat.Length] = -1;
ARowCopy[j][_cMat.Length + 1] = 1;
}
else
{
for (int k = 0; k < _cMat.Length; k++)
{
ARowCopy[j][k] = _aMat[j][k];
}
}
}
_aMat = ARowCopy;
//更新C值
double[] CRowCopy = new double[_cMat.Length + 2];
for (int j = 0; j < CRowCopy.Length; j++)
{
if (j < _cMat.Length)
{
CRowCopy[j] = _cMat[j];
}
else if (j == _cMat.Length)
{
CRowCopy[j] = 0;
}
else
{
CRowCopy[j] = M;
}
}
_cMat = CRowCopy;
//存储基变量下标
XBIndex[index++] = _cMat.Length - 1;
}
}
}
/// <summary>
/// 输出模型参数系数
/// </summary>
/// <param name="_cMat">目标函数系数向量</param>
/// <param name="_aMat">约束条件系数矩阵</param>
/// <param name="_bMat">约束条件资源向量</param>
private void OutParamater(double[] _cMat, double[][] _aMat, double[] _bMat, int[] _XBIndex)
{
Console.WriteLine("目标函数系数向量:");
for (int i = 0; i < _cMat.Length; i++)
{
Console.Write(_cMat[i] + " ");
}
Console.WriteLine();
Console.WriteLine("约束条件系数矩阵:");
for (int i = 0; i < _aMat.Length; i++)
{
for (int j = 0; j < _aMat[i].Length; j++)
{
Console.Write(AMat[i][j] + " ");
}
Console.Write(" " + _bMat[i]);
Console.WriteLine();
}
Console.WriteLine("基向量索引下标:");
for (int i = 0; i < _XBIndex.Length; i++)
{
Console.Write(_XBIndex[i] + " ");
}
Console.WriteLine();
}
}
}
第二部分
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
namespace SimplexMethod
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
private void button1_Click(object sender, EventArgs e)
{
Console.WriteLine("窗体调用成功!");
double[] c = new double[] { -2, -1 };
double[][] a = new double[][] { new double[] { 0, 5 }, new double[] { 6, 2 }, new double[] { 1, 1 } };
double[] b = new double[] { 15, 24, 5 };
int[] d = new int[] { -1, -1, -1 };
//double[] c = new double[] { -2,-3 };
//double[][] a = new double[][] { new double[] { -1, 1 }, new double[] { 1, 2 }, new double[] { 3, 1 } };
//double[] b = new double[] { 2,10, 15 };
//int[] d = new int[] { -1, -1,-1 };
//double[] c = new double[] { -4, -1 };
//double[][] a = new double[][] { new double[] { -1, 2 }, new double[] { 2, 3 }, new double[] { 1, -1 } };
//double[] b = new double[] { 4, 12, 3 };
//int[] d = new int[] { -1, -1, -1 };
//double[] c = new double[] { -6, -4 };
//double[][] a = new double[][] { new double[] { 1, 0 }, new double[] { 0, 1 }, new double[] { 3, 2 } };
//double[] b = new double[] { 3, 4, 11 };
//int[] d = new int[] { -1, -1, -1 };
SimplexMethodAlgorithm simplexMethodAlgorithm = new SimplexMethodAlgorithm(c, a, b,d);
simplexMethodAlgorithm.ExecuteAlgorithm();
List<double> solutionList = simplexMethodAlgorithm.solutionSet;
}
}
}