Двойственная задача/симплекс-метод [c] - C (СИ)
Формулировка задачи:
Дана была прямая задача, необходимо найти решение, а после по полученному решению найти решение двойственной ей задачи.
при данной системе ограничений:
Решить прямую задачу симплекс-методом получилось, вывод двойственной задачи, задачи в канонической форме тоже, но что делать дальше я не знаю, сказано, что решение двойственной задачи может быть найдено из формулы:
где D - матрица, составленная из векторов, входящих в последний базис, а
- вектор базисных переменных(последний получившийся в ходе преобразований)
Проблема-то, собственно, в том, что я не знаю, как получить обратную матрицу D, кроме того верно ее составив. И не помню, как верно перемножить вектор на матрицу.
И код не является универсальным, написан для конкретно этой задачи, и я не знаю, как теперь все пристроить.
В общем я совсем запутался, вполне возможно, забыл что-то ещё сделать или сделал не так.
#include <stdio.h> #include <stdlib.h> #include <float.h> #define WIDTH 5 #define HIGH 3 void printCanonic(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2); void printTask(const int arr1[], const double mat[][WIDTH], const double vect[], const int len); // вывод задачи void inverseMatrix(double **, double **, int); void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2); // вывод двойствеоон задачи double simplex(double A[HIGH][WIDTH], double A0[]); // 3-5 столбец - базис, 0 - вектор свободных членов void printSimplex(int [], int[], int [], int[], double[], double[], double[][HIGH], int, double); //вывод симлекс-таблицы int nminpos(double [], int); // минимальный отрицательный, позиция, -1 в случае отсутвия отриц. int pminpos(double [], double[], int); // позиция минимального неотрицательго частного double vectorMultiple(int [], double [], int); //скалярное void swap(int *, int *); const int f[5] = { 2, 3 }; // последние 3 для введенных переменных int main() { double A[HIGH][WIDTH] = { {-3,2,1,0,0}, { 1, -4, 0, 1, 0}, { 1, -1, 0, 0, 1} }; double b[HIGH] = { 6, 2, 5 }; printTask(f, A, b, 2); printCanonic(f,A,b,2,HIGH); printDuality(f,A,b,2,HIGH); simplex(A, b); return 0; } void printCanonic(const int arr1[], const double mat[][5], const double vect[], const int len1, const int len2) { int i, j; printf("##################################\n"); printf("# Task in canonical form #\n"); printf("##################################\n\n"); printf("f(x) = "); for (i = 0; i < len1; i++) printf("%c%dx%d", (i>0 && arr1[i]>0) ? '+' : ' ' ,arr1[i], i+1); printf("-> max\n"); printf("System of restrictions: \n"); for (i = 0; i < len2; i++) { for (j = 0; j < 5; j++) printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1); printf(" = %.lf\n", vect[i]); } putchar('\n'); } void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2) { int i, j; printf("##################################\n"); printf("# Twofold task #\n"); printf("##################################\n\n"); printf("g(y) = "); for (i = 0; i < len2; i++) printf("%c%.lfy%d", (i>0 && vect[i]>0) ? '+' : ' ' ,vect[i], i+1); printf("-> min\n"); printf("System of restrictions: \n"); for (i = 0; i < len1; i++) { for (j = 0; j < len2; j++) printf("%c%.lfx%d", (j>0 && mat[j][i]>=0) ? '+' : ' ' ,mat[j][i], j+1); printf(" >= %d\n", arr1[i]); } putchar('\n'); } void printTask(const int arr[], const double mat[][WIDTH], const double vect[], const int len1) { int i, j; printf("##################################\n"); printf("# Straight task #\n"); printf("##################################\n\n"); printf("f(x) = "); for (i = 0; i < len1; i++) printf("%c%dx%d", (i>0 && arr[i]>0 && arr[i] != 1) ? '+' : ' ' ,arr[i], i+1); printf("-> max\n"); printf("System of restrictions: \n"); for (i = 0; i < HIGH; i++) { for (j = 0; j < len1; j++) printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1); printf(" <= %.lf\n", vect[i]); } putchar('\n'); } double simplex(double A[HIGH][WIDTH], double A0[]) { int Cbasis[HIGH] = { 0, 0, 0 }; int basis[HIGH] = {3, 4, 5}; int notbase[2] = {1,2}; int Cj[2] = { f[0], f[1] }; double fstr[2]/*, x[WIDTH] = { 0, 0, A0[0], A0[1], A0[2]}*/; double matr[WIDTH][HIGH], commonmatr[WIDTH][HIGH], Q, **D, **D1; int r, s, i, j, counter = 0; for(i=0; i<WIDTH; i++) for (j=0; j<HIGH; j++) matr[i][j] = A[j][i]; for (i=0; i<2;i++) fstr[i] = vectorMultiple(Cbasis, matr[i], HIGH) - f[i]; Q = vectorMultiple(Cbasis, A0, HIGH); printf("\nSimplex-matrix:\n"); printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q); while ((r = nminpos(fstr, 2)) != -1) { counter++; if ((s = pminpos(matr[r], A0, HIGH)) == -1) return -1; for(i=0; i<WIDTH; i++) for (j=0; j<HIGH; j++) commonmatr[i][j] = matr[i][j]; matr[r][s] = 1/commonmatr[r][s]; for (i = 0; i < HIGH; i++) if (i!=s) matr[r][i] = -commonmatr[r][i] / commonmatr[r][s]; for (i = 0; i < 2; i++) if (i!=r) matr[i][s] = commonmatr[i][s] / commonmatr[r][s]; for (i = 0; i < HIGH; i++) for (j=0; j<2; j++) if ((i != r) && (j!=s)) matr[i][j] = (commonmatr[i][j] * commonmatr[r][s] - commonmatr[r][j] * commonmatr[i][s])/commonmatr[r][s]; for (i = 0; i < HIGH; i++) { if (i!=s) A0[i] = (A0[i] * commonmatr[r][s] - commonmatr[r][i] * A0[s])/commonmatr[r][s]; } A0[s] /= commonmatr[r][s]; //f-строка for (i = 0; i < 2; i++) { if (i!=r) fstr[i] = (fstr[i] * commonmatr[r][s] - commonmatr[i][s] * f[r])/commonmatr[r][s]; } fstr[r] /= (-commonmatr[r][s]); swap(&Cbasis[s], &Cj[r]); swap(&basis[s], ¬base[r]); Q = vectorMultiple(Cbasis, A0, HIGH); printf("\nSimplex-matrix after %d iteration:\n", counter); printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q); } return Q; } void swap(int *x, int *y) { int temp = *x; *x = *y; *y = temp; } int pminpos(double A[], double A0[], int length) { int i, minp = -1; double min; min = DBL_MAX; for (i = 0; i < length; i++) { if ((A[i] > 0) && (A0[i]/A[i] < min)) { min = A0[i]/A[i]; minp = i; } } return minp; } int nminpos(double A[], int length) { int minp = -1, i; double min = 0; for (i = 0; i < length; i++) { if (A[i] < min) { min = A[i]; minp = i; } } return minp; } void printSimplex(int Cbas[], int bas[], int cj[], int notb[], double fs[], double a0[], double matrix[][HIGH], int clmn, double Q) { int i; printf("_____|__Cj__|__%-5d|__%-5d|_____\n", cj[0], cj[1]); printf("Basis|______|_x%-5d|_x%-5d|__A0_\n", notb[0], notb[1]); for (i = 0; i < HIGH; i++) printf("__%-3d|_x%-4d|%7.2lf|%7.2lf|%5.2lf\n", Cbas[i], bas[i], matrix[0][i], matrix[1][i], a0[i]); printf("_____|___f__|%7.2lf|%7.2lf|%5.2lf\n", fs[0], fs[1], Q); printf("_____|______|_DELTA1|_DELTA2|__Q__\n"); } double vectorMultiple(int a1[], double a2[], int length) { int i; double sum = 0; for (i = 0; i < length; i++) sum += a1[i]*a2[i]; return sum; }
Решение задачи: «Двойственная задача/симплекс-метод [c]»
textual
Листинг программы
#include <stdio.h> #include <stdlib.h> #include <float.h> #define WIDTH 5 #define HIGH 3 void printCanonic(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2); void printTask(const int arr1[], const double mat[][WIDTH], const double vect[], const int len); void inverseMatrix(double **, int); void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2); double simplex(double A[HIGH][WIDTH], double A0[]); // 3-5 базис, 0 - вектор свободных членов void printSimplex(int [], int[], int [], int[], double[], double[], double[][HIGH], int, double); int nminpos(double [], int); int pminpos(double [], double[], int); double vectorMultiple(int [], double [], int); void swap(int *, int *); void matrXvect(double **, double *, int); const int f[5] = { 2, 3 }; int main() { double A[HIGH][WIDTH] = { {-3,2,1,0,0}, { 1, -4, 0, 1, 0}, { 1, -1, 0, 0, 1} }; double b[HIGH] = { 6, 2, 5 }, result; printTask(f, A, b, 2); printCanonic(f,A,b,2,HIGH); printDuality(f,A,b,2,HIGH); result = simplex(A, b); if (result!=-1) printf("\nMaximim is %.2lf\nResult is correct.\n", result); else printf("\nInclorrect result/undecidability problem\n"); return 0; } void printCanonic(const int arr1[], const double mat[][5], const double vect[], const int len1, const int len2) { int i, j; printf("##################################\n"); printf("# Task in canonical form #\n"); printf("##################################\n\n"); printf("f(x) = "); for (i = 0; i < len1; i++) printf("%c%dx%d", (i>0 && arr1[i]>0) ? '+' : ' ' ,arr1[i], i+1); printf("-> max\n"); printf("System of restrictions: \n"); for (i = 0; i < len2; i++) { for (j = 0; j < 5; j++) printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1); printf(" = %.lf\n", vect[i]); } putchar('\n'); } void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2) { int i, j; printf("##################################\n"); printf("# Twofold task #\n"); printf("##################################\n\n"); printf("g(y) = "); for (i = 0; i < len2; i++) printf("%c%.lfy%d", (i>0 && vect[i]>0) ? '+' : ' ' ,vect[i], i+1); printf("-> min\n"); printf("System of restrictions: \n"); for (i = 0; i < len1; i++) { for (j = 0; j < len2; j++) printf("%c%.lfx%d", (j>0 && mat[j][i]>=0) ? '+' : ' ' ,mat[j][i], j+1); printf(" >= %d\n", arr1[i]); } putchar('\n'); } void printTask(const int arr[], const double mat[][WIDTH], const double vect[], const int len1) { int i, j; printf("##################################\n"); printf("# Straight task #\n"); printf("##################################\n\n"); printf("f(x) = "); for (i = 0; i < len1; i++) printf("%c%dx%d", (i>0 && arr[i]>0 && arr[i] != 1) ? '+' : ' ' ,arr[i], i+1); printf("-> max\n"); printf("System of restrictions: \n"); for (i = 0; i < HIGH; i++) { for (j = 0; j < len1; j++) printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1); printf(" <= %.lf\n", vect[i]); } putchar('\n'); } double simplex(double A[HIGH][WIDTH], double A0[]) { int Cbasis[HIGH] = { 0, 0, 0 }; int basis[HIGH] = {3, 4, 5}; int notbase[2] = {1,2}; int Cj[2] = { f[0], f[1] }; double fstr[2]/*, x[WIDTH] = { 0, 0, A0[0], A0[1], A0[2]}*/; double matr[WIDTH][HIGH], commonmatr[WIDTH][HIGH], Q, Qy, D[3][3], *x[3], basdouble[HIGH], commonA0[HIGH]; int r, s, i, j, counter = 0; for (j=0; j<HIGH; j++) commonA0[j] = A0[j]; for(i=0; i<WIDTH; i++) for (j=0; j<HIGH; j++) matr[i][j] = A[j][i]; for (i=0; i<2;i++) fstr[i] = vectorMultiple(Cbasis, matr[i], HIGH) - f[i]; Q = vectorMultiple(Cbasis, A0, HIGH); printf("\nSimplex-matrix:\n"); printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q); while ((r = nminpos(fstr, 2)) != -1) { counter++; if ((s = pminpos(matr[r], A0, HIGH)) == -1) return -1; for(i=0; i<WIDTH; i++) for (j=0; j<HIGH; j++) commonmatr[i][j] = matr[i][j]; matr[r][s] = 1/commonmatr[r][s]; for (i = 0; i < HIGH; i++) if (i!=s) matr[r][i] = -commonmatr[r][i] / commonmatr[r][s]; for (i = 0; i < 2; i++) if (i!=r) matr[i][s] = commonmatr[i][s] / commonmatr[r][s]; for (i = 0; i < HIGH; i++) for (j=0; j<2; j++) if ((i != r) && (j!=s)) matr[i][j] = (commonmatr[i][j] * commonmatr[r][s] - commonmatr[r][j] * commonmatr[i][s])/commonmatr[r][s]; for (i = 0; i < HIGH; i++) { if (i!=s) A0[i] = (A0[i] * commonmatr[r][s] - commonmatr[r][i] * A0[s])/commonmatr[r][s]; } A0[s] /= commonmatr[r][s]; //f-строка for (i = 0; i < 2; i++) { if (i!=r) fstr[i] = (fstr[i] * commonmatr[r][s] - commonmatr[i][s] * f[r])/commonmatr[r][s]; } fstr[r] /= (-commonmatr[r][s]); swap(&Cbasis[s], &Cj[r]); swap(&basis[s], ¬base[r]); Q = vectorMultiple(Cbasis, A0, HIGH); printf("\nSimplex-matrix after %d iteration:\n", counter); printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q); } for (i=0; i<3; i++) for (j=0;j<3;j++) D[i][j] = A[j][basis[i]-1]; printf("\nMatrix D of last basis columns:\n"); for (i=0; i<3; i++) { x[i] = D[i]; for (j=0;j<3;j++) printf("%.2lf ", x[i][j] = D[i][j]); putchar('\n');} inverseMatrix(x, 3); printf("\nInvert matrix D^-1:\n"); for (i=0; i<3; i++) { for (j=0;j<3;j++) printf("%.2lf ", x[i][j]); putchar('\n');} for (i=0; i<3; i++) basdouble[i] = Cbasis[i]; matrXvect(x, basdouble, 3); printf("\nVector y = ("); for (i=0; i<3; i++) printf("%.2lf, ",basdouble[i]); printf(")\n"); printf("\nValue of duality funclion:\n"); printf("g(y) = "); Qy = 0; for (i = 0; i < HIGH; i++) { Qy+=commonA0[i]*basdouble[i]; printf("%c%.lfy%d", (i>0 && commonA0[i]>0) ? '+' : ' ' ,commonA0[i], i+1); } printf(" = %.2lf\n", Qy); if (Q!=Qy) return -1; return Q; } void swap(int *x, int *y) { int temp = *x; *x = *y; *y = temp; } int pminpos(double A[], double A0[], int length) { int i, minp = -1; double min; min = DBL_MAX; for (i = 0; i < length; i++) { if ((A[i] > 0) && (A0[i]/A[i] < min)) { min = A0[i]/A[i]; minp = i; } } return minp; } int nminpos(double A[], int length) { int minp = -1, i; double min = 0; for (i = 0; i < length; i++) { if (A[i] < min) { min = A[i]; minp = i; } } return minp; } void printSimplex(int Cbas[], int bas[], int cj[], int notb[], double fs[], double a0[], double matrix[][HIGH], int clmn, double Q) { int i; printf("_____|__Cj__|__%-5d|__%-5d|_____\n", cj[0], cj[1]); printf("Basis|<br>|_x%-5d|_x%-5d|__A0_\n", notb[0], notb[1]); for (i = 0; i < HIGH; i++) printf("__%-3d|_x%-4d|%7.2lf|%7.2lf|%5.2lf\n", Cbas[i], bas[i], matrix[0][i], matrix[1][i], a0[i]); printf("_____|___f__|%7.2lf|%7.2lf|%5.2lf\n", fs[0], fs[1], Q); printf("_____|<br>|_DELTA1|_DELTA2|__Q__\n"); } double vectorMultiple(int a1[], double a2[], int length) { int i; double sum = 0; for (i = 0; i < length; i++) sum += a1[i]*a2[i]; return sum; } void inverseMatrix(double **A, int n) { ///Метод гауссса-йордана int i,j,k,m; double akk, aik, r, E[HIGH][HIGH] = { {1,0,0}, {0,1,0}, {0,0,1} }; m = 0; for (k=0; k <= n-1; k++) { if ((A[k][k]) == 0) { m = k +1; while (m<=n && A[m][k]==0) m++; if (m > n)// { printf("Matrix is degenerate\n"); //return -1 } for (i = 0; i<n; i++) { r = A[m][i]; A[m][i] = A[k][i]; A[k][i] = r; r = E[m][i]; E[m][i] = E[k][i]; E[k][i] = r; } } akk = A[k][k]; //A[k][k]/=akk; E[k][k]/=akk; for (i = 0; i < n; i++) if (i!=k) { aik = A[i][k]/akk; for (j = k; j < n; j++) { A[i][j] = A[i][j] - aik*A[k][j]; E[i][j] = (E[i][j] - aik*E[k][j])/A[i][i]; } }; }; for (i = 0; i < n; i++) for (j = 0; j < n; j++) A[i][j] = E[i][j]; } void matrXvect(double **A, double *c, int n) { double ii; int i, j; for (i = 0; i < n; i++) { ii = 0; for (j = 0; j < n; j++) ii += A[j][i]*c[i]; c[i] = ii; } }
ИИ поможет Вам:
- решить любую задачу по программированию
- объяснить код
- расставить комментарии в коде
- и т.д