//файл включения для векторных и матричных объектов
#include "matrix.h"
/*определяем типы "действительная матрица" и "действительный вектор"*/
typedef matrix<double> dmatrix;
typedef vector<double> dvector;
/*процедура перестановки векторов матрицы для выбора главного элемента*/
void rotate(dmatrix &x, long index)
{
/*нам необходимо, чтобы в процедуре прямого хода на главной диагонали не было нулевых составляющих; если же таковые есть, находим "ниже" по матрице вектор-уравнение с ненулевой соответствующей компонентой и ставим его на место текущего*/
for(long i=index;i<x.getm();i++)
if(x[i][index])
{
dvector temp=x[i];
x[i]=x[index];
x[index]=temp;
}
}
/*функция, выполняющая прямой ход Гаусса с попутным подсчётом количества линейно-независимых векторов данной матрицы (ранга матрицы)*/
long rank(dmatrix &x)
{
long i;
//перебираем вектора матрицы
for(i=0;i<x.getm();i++)
{
if(!x[i][i])/*если диагональный элемент равен нулю*/
rotate(x,i);//меняем строки матрицы
/*делаем нулевым столбец под текущим диагональным элементом*/
for(long j=i+1;j<x.getm();j++)
if(x[i][i])
x[j]+=(-x[j][i]/x[i][i])*x[i];
}
//создаём пустой вектор заданной размерности
dvector empty=x.getn();
/*количество нулевых векторов в матрице после прямого хода и определяет ранг данной матрицы*/
for(i=0;i<x.getm();i++)
if(x[i]==empty)/*сравниваем текущий вектор с нулевым*/
break;
return i;//возвращаем ранг
}
/*процедура обратного хода Гаусса без нормирования элементов главной диагонали*/
void diagonalize(dmatrix &x)
{
for(long i=x.getm()-1;i>=0;i--)
for(long j=i-1;j>=0;j--)
if(x[i][i])
x[j]+=(-x[j][i]/x[i][i])*x[i];
}
/*Процедура обмена двух столбцов матрицы - перестановки местами переменных с попутным переразрешением полученной системы. Параметрами её являются матрица ОЗЛП, вектор линейной формы, номера обмениваемых переменных и вектор, в котором учитывается производимая перестановка*/
void swap(dmatrix &x, dvector &l, long sw1, long sw2, dvector &swapvalues)
{
//временный вектор для обмена значениями
dvector save(x.getm()+2);
//сохраняем первый столбец
for(long i=0;i<x.getm();i++)
save[i]=x[i][sw1];
/*сохраняем элемент линейной формы, соответствующий первому столбцу*/
save[x.getm()]=l[sw1];
//сохраняем номер переставляемого столбца и т.д.
save[x.getm()+1]=swapvalues[sw1];
for(long i=0;i<x.getm();i++)
x[i][sw1]=x[i][sw2];
l[sw1]=l[sw2];
swapvalues[sw1]=swapvalues[sw2];
for(long i=0;i<x.getm();i++)
x[i][sw2]=save[i];
l[sw2]=save[x.getm()];
swapvalues[sw2]=save[x.getm()+1];
rank(x);//осуществляем прямой и
diagonalize(x);//обратный ход Гаусса
for(long i=0;i<x.getm();i++)
x[i]*=(1/x[i][i]);
/*подставляем значений переменных в линейную форму, определяя её коэффициенты в данной перестановке*/
for(long i=0;i<x.getm();i++)
l+=-l[i]*x[i];
}
/*процедура определения опорного решения принимает три параметра - матрицу ОЗЛП, вектор линейной формы и вектор перестановок*/
/*ищем уравнение с отрицательным свободным членом (заметим, что в нашей записи все переменные находятся в левой части уравнения, поэтому коэффициенты при них следует брать с обратным знаком)*/
for(i=0;i<x.getm();i++)
if(-x[i][x.getn()-1]<0)
break;
if(i==x.getm())
return;/*опорное решение найдено - все свободные члены неотрицательны */
/*в строке, где мы нашли отрицательный свободный член, среди коэффициентов пытаемся найти положительный*/
for(j=x.getm();j<x.getn()-1;j++)
if(x[i][j]<0)
break;
/*если положительный коэффициент мы найти не смогли, то нарушается принцип неотрицательности переменных, что недопустимо*/
if(j==x.getn()-1)
throw xmsg("Допустимых решений нет");
/*среди строк, в которых в определённом (разрешающем) столбце отрицательному свободному члену соответствует положительный коэффициент, выбираем ту, в которой отношение -x[k][n-1]/x[k][j] принимает наименьшее значение*/
for(k=0;k<x.getm();k++)
if(sign(x[k][j])==sign(-x[k][x.getn()-1])
&&x[k][j])
{
s=k;
otn=-x[k][x.getn()-1]/x[k][j];
break;
}
for(;k<x.getm();k++)
if(sign(x[k][j])==sign(-x[k][x.getn()-1])
&&x[k][j])
if(otn>-x[k][x.getn()-1]/x[k][j])
s=k,
otn=-x[k][x.getn()-1]/x[k][j];
/*j - разрешающий столбец, s - разрешающая строка. Меняем в нашей матрице столбцы с соответствующими индексами */
swap(x, l, s, j, swapvalues);
//снова пытаемся определить опорное решение
getoporsol(x,l,swapvalues);
}
/*если вектор оптимального решения существует, то данная функция его определит, используя матрицу ОЗЛП, линейную форму и вектор перестановок*/
/*если таковых нет, полученное решение является оптимальным*/
if(i==l.getm()-1)
{
//формируем вектор результата
dvector res=l.getm()-1;
/*используя массив перестановок, записываем значения компонент вектора решения на их прежние позиции*/
for(long i=0;i<x.getm();i++)
res[swapvalues[i]]=-x[i][x.getn()-1];
return res;//возвращаем оптимальное решение
}
/*иначе - определяем номер "опасного" уравнения, в котором увеличение переменной, соответствующей отрицательному коэффициенту линейной формы, может привести к нарушению условия неотрицательности переменных */
for(j=0;j<x.getm();j++)
if(-x[j][i]<0)
break;
/*если такого уравнения нет, то данную переменную можно увеличивать безгранично, а, значит, оптимального решения ОЗЛП не существует*/
if(j==x.getm())
throw xmsg("Линейная форма не ограничена снизу");
/*если же такие уравнения нашлись, определяем среди них ту переменную, которая при увеличении данной наиболее быстро обратится в нуль*/
for(k=0;k<x.getm();k++)
if(-x[k][i]<0)//опасное уравнение
{
s=k;
otn=-x[k][x.getn()-1]/x[k][i];
break;
}
for(;k<x.getm();k++)
if(-x[k][i]<0)//опасное уравнение
if(otn>-x[k][x.getn()-1]/x[k][i])
s=k,//наиболее угрожаемая переменная
otn=-x[k][x.getn()-1]/x[k][i];
/*меняем местами переменную при отрицательном коэффициенте линейной формы с наиболее угрожаемой */
swap(x,l,s,i,swapvalues);
/*после замены вновь пытаемся найти оптимальное решение*/
return getoptsol(x,l,swapvalues);
}
/*функция для решения ОЗЛП, заданной основной матрицей (левой частью), столбцом свободных членов правой части и линейной формой без св. члена */
dvector simplex(dmatrix a, dvector b, dvector c)
{
if(a.getm()!=b.getm())
throw xmsg("Количество уравнений должно совпадать");
if(a.getn()!=c.getm())
throw xmsg("Количество неизвестных должно совпадать");
/*матрица для хранения системы уравнений, заданной левой и правой частями*/
dmatrix ab(a.getm(),a.getn()+1);
dvector swapvalues=c.getm();/*вектор учёта обмена переменных*/
//вначале порядок переменных не нарушен
for(long i=0;i<swapvalues.getm();i++)
swapvalues[i]=i;
//формируем матрицу СЛАУ из левой и правой частей
for(long i=0,j;i<a.getm();i++)
{
for(j=0;j<a.getn();j++)
ab[i][j]=a[i][j];
ab[i][j]=-b[i];
}
long r;
//определяем ранги основной и расширенной матрицы
if((r=rank(a))!=rank(ab))
throw xmsg("Система несовместна - ранги основной и расширенной матриц не совпадают");
/*в случае совместной системы ранги этих матриц совпадают, поэтому мы можем отбросить уравнения-следствия*/
dmatrix x(r,ab.getn());
/*в новую матрицу ОЗЛП переписываем только линейно-независимые уравнения*/
for(long i=0;i<r;i++)
x[i]=ab[i];
/*добавляем в вектор линейной формы пока нулевой свободный член*/
dvector l(c.getm()+1);
for(long i=0;i<c.getm();i++)
l[i]=c[i];
//выполняем обратный ход Гаусса
diagonalize(x);
for(long i=0;i<x.getm();i++)
x[i]*=(1/x[i][i]);
/*подстановка значений базисных переменных в линейную форму*/
for(long i=0;i<x.getm();i++)
l+=-l[i]*x[i];
//поиск опорного решения
getoporsol(x,l,swapvalues);
//возвращаем оптимальное решение
return getoptsol(x,l,swapvalues);
}
//тестовая программа работает с задачей вида:
//-5*x1-x2+2*x3<=2
//-x1+x3+x4<=5
//-3*x1+5*x4<=7 L=5*x1-2*x3
//x1>=0
//x2>=0
//x3>=0
//x4>=0
//==>
//5*x1+x2-2*x3>=-2
//x1-x3-x4>=-5
//3*x1-5*x4>=-7 L=5*x1-2*x3
//x1>=0
//x2>=0
//x3>=0
//x4>=0
//==>
//5*x1+x2-2*x3+2>=0
//x1-x3-x4+5>=0
//3*x1-5*x4+7>=0 L=5*x1-2*x3
//x1>=0
//x2>=0
//x3>=0
//x4>=0
//==>
//y1=5*x1+x2-2*x3+2
//y2=x1-x3-x4+5
//y3=3*x1-5*x4+7 L=5*x1-2*x3
//y1>=0
//y2>=0
//y3>=0
//x1>=0
//x2>=0
//x3>=0
//x4>=0
//==>
//-y1+5*x1+x2-2*x3+2=0
//-y2+x1-x3-x4+5=0
//-y3+3*x1-5*x4+7=0 L=5*x1-2*x3
//y1>=0
//y2>=0
//y3>=0
//x1>=0
//x2>=0
//x3>=0
//x4>=0
//
void main()
{
double mtr[]=
{
-1,0,0, 5,1,-2,0,
0,-1,0, 1,0,-1,-1,
0,0,-1, 3,0,0,-5
};
double vec1[]={-2, -5, -7};
double vec2[]={0,0,0,5,0,-2,0};
dmatrix a(3,7,mtr), b(3,vec1), c(7,vec2);
dvector res=simplex(a,b,c);
cout<<res<<endl;
/*значение линейной формы - скалярное произведение вектора переменных на вектор коэффициентов, т.е. сумма соответствующих произведений*/
