Код:
#include <vector>
#include <string>
#include <math.h>
using namespace std;
#define SIMPLEX_DONE 0 // оптимизация успешно завершена
#define SIMPLEX_NO_SOLUTION 1 // задача не имеет решения (не удается найти базис)
#define SIMPLEX_NO_BOTTOM 2 // решения нет, т.к. линейная форма не ограничена снизу
#define SIMPLEX_NEXT_STEP 3 // для получения решения нужно сделать еще хотя бы один шаг
#define MAX_VAL 0.1e-12 //точность (значение, удовлетворяющее -MAX_VAL < X < MAX_VAL считается нулем)
#define MaxExtended 1.1e+100 //Максимальное значение
//1.7E+308
//extended -3.6 x 10^4951 до 1.1 x 10^4932
//MaxExtended 1.1 E+ 4932
typedef enum { Equal,Less,Greater } TOperation;
typedef struct TConstrain
{
vector<double> A;
double B;
TOperation Sign;
bool isT;
} TConstrain;
class TSimplex
{
public:
int M, N; // M - число строк, N - число столбцов
int RealN; // реальное число переменных, изначально вошедших в задачу
vector<TConstrain> Cons;
vector<double> C;
double L;
vector<int> Basis;
bool Max;
TSimplex::TSimplex()
{
}
//Конструктор
TSimplex::TSimplex(vector<double> _C, bool Maximize)
{
int j;
N = _C.size();
RealN = N;
M = 0;
C.resize(N);
Max = Maximize;
if ((!Maximize))
for (j = 0; j < N ; j++)
C[j] = -_C[j];
else
for (j = 0; j < N ; j++)
C[j] = _C[j];
Max = Maximize;
L = 0;
}
//Если num лежит внутри области точности значений, возвращает 0
double TSimplex::DoPrec(double num)
{
if ((num < MAX_VAL) && (num > -MAX_VAL))
num = 0;
return num;
}
//Добавление ограничения вида Ai[1]*x1+Ai[2]*x2+....+Ai[n]*xn <=> Bi
void TSimplex::AddCons(double _B, vector<double> _A, TOperation Sign)
{
int j;
if ((_A.size() > N))
SetAllLengths(_A.size());
M++;
Cons.resize(M);
Cons[M - 1].B = _B;
Cons[M - 1].Sign = Sign;
Cons[M - 1].A.resize(N);
for (j = 0; j < _A.size() ; j++)
Cons[M - 1].A[j] = _A[j];
if (_A.size() < N)
for (j = _A.size(); j < N ; j++)
Cons[M-1].A[j] = 0;
}
//Суммирование строки 1 со строкой 2, домноженной на коэффициент Value
void TSimplex::AddString(int Num1,int Num2,double Value)
{
int j;
for (j = 0; j <= N - 1; j++)
Cons[Num1].A[j] = Cons[Num1].A[j] + Cons[Num2].A[j] * Value;
Cons[Num1].B = Cons[Num1].B + Cons[Num2].B * Value;
}
//
bool TSimplex::CheckBasis(void)
{
int i, j, k;
bool f;
Basis.resize(M);
for (i = 0; i <M; i++)
Basis[i] = -1;
for (j = 0; j <= N - 1; j++)
{
f = true;
k = -1;
i = 0;
while ((f && (i < M)))
{
if (((Cons[i].A[j] != 0) && (Cons[i].A[j] != 1)))
f = false;
if ((Cons[i].A[j] == 1))
{
if ((k == -1)) k = i;
else f = false;
}
i+=1;
}
if ((f && (k != -1)))
Basis[k] = j;
}
f = true;
for (i = 0; i <= M - 1; i++)
f = f && (Basis[i] != -1);
return f;
}
void TSimplex::CreateBasis(TSimplex Simplex)
{
int i, j;
M = Simplex.M;
N = Simplex.N;
RealN = Simplex.RealN;
L = 0;
Cons.resize(M);
Basis.resize(M);
C.resize(N);
for (i = 0; i < N ; i++)
C[i] = 0;
for (i = 0; i < M ; i++)
{
Cons[i].A.resize(N);
for (j = 0; j < N; j++)
Cons[i].A[j] = Simplex.Cons[i].A[j];
Cons[i-1].B = Simplex.Cons[i-1].B;
Cons[i-1].Sign = Equal;
Cons[i-1].isT = false;
}
for (i = 0; i < M; i++)
{
if ((Simplex.Basis[i] != -1))
Basis[i] = Simplex.Basis[i];
else
{
SetAllLengths(N + 1);
for (j = 0; j <= M - 1; j++)
Cons[j].A[N - 1] = 0;
Cons[i].A[N - 1] = 1;
Cons[i].isT = true;
C[N - 1] = 0;
for (j = 0; j <= Simplex.N - 1; j++)
C[j] = C[j] + Simplex.Cons[i].A[j];
L = L + Cons[i].B;
}
}
}
void TSimplex::Free()
{
C.resize(0);
Basis.resize(0);
Cons.resize(0);
M = 0;
N = 0;
RealN = 0;
}
double TSimplex::GetMin()
{
//int i;
if ((Max))
return -L;
else
return L;
}
vector<double> TSimplex::GetSolution()
{
vector<double> Solution;
int i, j;
Solution.resize(RealN);
for (j = 0; j < RealN; j++)
{
Solution[j] = 0;
i = 0;
while (((i < M) && (Basis[i] != j)))
i++;
if (((Basis[i] == j) && (i < M)))
Solution[j] = Cons[i].B;
}
return Solution;
}
void TSimplex::MulString(int Number,double Value)
{
int j;
for (j = 0; j < N; j++)
Cons[Number].A[j] = Cons[Number].A[j] * Value;
Cons[Number].B = Cons[Number].B * Value;
}