' . "\n"; $str .= '' . "\n"; $str .= '' . "\n"; /* Build head */ $str .= "\n"; $str .= '\n"; $str .= "Simplex method\n"; $str .= '' . "\n"; $str .= '' . "\n"; $str .= '' . "\n"; $str .= '' . "\n"; $str .= '' . "\n"; $str .= '' . "\n"; $str .= '' . "\n"; $str .= '' . "\n"; $str .= ''; $str .= "\n\n"; /* Build body */ $str .= '

Simplex Method

' . "\n"; $str .= '
' . "\n"; $str .= '
Select the data file from your computer (separator semicolon, last row objective function, first column right-hand sides: cf. example):' . "
\n"; $str .= '' . "\n"; $output = isset($_REQUEST['output']) || !isset($_REQUEST['submit']); $str .= ($output) ? '' . "\n" : '' . "\n"; if (isset($_REQUEST['prog'])) $prog = $_POST['prog']; else $prog = 1; $str .= ($prog == 1) ? '' . "\n" : '' . "\n"; $str .= ($prog == 2) ? '' . "\n" : '' . "\n"; if (isset($_REQUEST['rule'])) $rule = $_POST['rule']; else $rule = 1; $str .= ($rule == 1) ? '' . "\n" : '' . "\n"; $str .= ($rule == 2) ? '' . "\n" : '' . "\n"; $str .= ($rule == 3) ? '' . "\n" : '' . "\n"; if (isset($_REQUEST['opt'])) $opt = $_POST['opt']; else $opt = 1; $str .= ($opt == 1) ? '' . "\n" : '' . "\n"; $str .= ($opt == 2) ? '' . "\n" : '' . "\n"; if (isset($_REQUEST['op'])) $op = $_POST['op']; else $op = 1; $str .= ($op == 1) ? '' . "\n" : '' . "\n"; $str .= ($op == 2) ? '' . "\n" : '' . "\n"; $red = isset($_REQUEST['red']) || !isset($_REQUEST['submit']); $str .= ($red) ? '' . "\n" : '' . "\n"; $nor1 = isset($_REQUEST['nor1']); $str .= ($nor1) ? '' . "\n" : '' . "\n"; $scal = isset($_REQUEST['scal']); $str .= ($scal) ? '' . "\n" : '' . "\n"; $nor2 = isset($_REQUEST['nor2']); $str .= ($nor2) ? '' . "\n" : '' . "\n"; $file = ($_FILES && $_FILES['userfile']['type'] == "text/plain") ? $_FILES['userfile']['name'] : ((isset($_REQUEST['file'])) ? $_POST['file'] : "bh_de_data.txt"); unset($_FILES); $str .= '' . "\n"; $str .= '
' . "\n"; /* Definitions and functions */ define('INFINITY', INF); define('EPSILON', 1E-10); function format($val) { /* Format single value */ if ($val == "-0") $val = 0; return (is_finite(floatval($val))) ? strtoupper(str_replace('.', ',', $val)) : (($val > 0) ? "∞" : "-∞"); } function swap(&$left, &$right) { /* Swap variables */ $tmp = $left; $left = $right; $right = $tmp; } function zero(&$a, &$i, &$j, $m, $n, &$solve) { /* Set small values to zero and check if there are infinite values */ for ($i = 0; $i < $m; $i++) for ($j = 0; $j < $n; $j++) if (!is_finite(floatval($a[$i][$j]))) return $solve = false; else if (abs($a[$i][$j]) < EPSILON) $a[$i][$j] = 0; return $solve = true; } function normalise(&$a, &$no, &$i, &$j, $m, $n) { /* Normalise constraints line by line */ for ($i = 1; $i < $m; $i++) { $sca = 0; for ($j = 1; $j < $n; $j++) $sca += $a[$i][$j] * $a[$i][$j]; $sca = Sqrt($sca); for ($j = 0; $j < $n; $j++) $a[$i][$j] /= $sca; $no[$i] = 1 / $sca; } } function scale(&$a, &$e, &$i, &$j, $m, $n) { /* Scale constraints column by column */ for ($j = 1; $j < $n; $j++) { $min = INFINITY; for ($i = 1; $i < $m; $i++) { $obj = abs($a[$i][$j]); if ($obj > 0 && $obj < $min) $min = $obj; } $min = 1 / $min; for ($i = 0; $i < $m; $i++) $a[$i][$j] *= $min; $e[$j] = $min; } } function redundant(&$a, &$bv, &$mbv, &$cb, &$cn, $m, $n) { /* Check whether constraint is <= 0 and push it then increase counter by one if it is new */ for ($i = 1; $i < $m; $i++) { $test = $bv[$i]; if ($mbv[$test]) { for ($j = 1; $j < $n && $a[$i][$j] <= 0; $j++); if ($j == $n) { $mbv[$test] = false; if ($test < $n) $cn++; else $cb++; } } } } function trans(&$a, $k, $l, $m, $n, $bv, $nbv, $opt, $prog, &$solve, &$str) { /* Boundary of computation was transcended */ $solve = false; tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); $str .= "
The boundary of computation was transcended.
\n"; } function check(&$a, &$m, &$n, &$bv, &$nbv, $opt, $prog, &$solve, &$str, $lines) { /* Check wrong characters */ $input = str_replace('
', '', $lines) . ";"; $rest = preg_filter('/[\d\+\-\.;eE\s]/', '', $input); if ($rest != '') { preg_match_all('/./', $rest, $signs); $error = $sign = ''; sort($signs[0]); for ($i = 0; $i < count($signs[0]); $i++) { if ($signs[0][$i] != $sign) { $sign = $signs[0][$i]; $error .= "$sign, "; } } if ($sign == $signs[0]) $str .= "The character $sign is not permitted.\n"; else $str .= "The following characters are not permitted: " . substr($error, 0, strlen($error) - 2) . ".\n"; return false; } /* Check numbers */ if (preg_match('/^((\+|-)?(\d{0,16}\.\d{1,16}|\d{1,16})((e|E)(\+|-)?\d+)?([;\n\r]|\r\n))+$/', $input)) { $input = preg_filter('/[\n\r]/', '$', $input); $input = strtoupper(str_replace('$$', '$', $input)); if (substr($input, strlen($input) - 1, 1) == "$") $input = substr($input, 0, strlen($input) - 1); $rex = '/(\-?\d{0,16}\.\d{1,16}(E\-?\d+)?|\-?\d+(E\d+)?(\/\d+(E\d+)?)?)/'; $rows = explode('$', $input); preg_match_all($rex, $rows[0], $hlp); $a[1] = $hlp[0]; $n = count($a[1]); $m = count($rows); /* Check matrix */ for ($k = 2; $k < $m; $k++) { preg_match_all($rex, $rows[$k - 1], $hlp); $a[$k] = $hlp[0]; $len = count($a[$k]); if ($n != $len) { if ($len == 1) $str .= "The matrix has in row $k unfortunately one column instead $n.\n"; else $str .= "The matrix has in row $k unfortunately $len columns instead $n.\n"; return false; } } if ($m > 1) { preg_match_all($rex, $rows[$k - 1], $hlp); $a[0] = $hlp[0]; $len = count($a[0]); if ($n != $len) { if ($len == 1) $str .= "The matrix has in row $k unfortunately one column instead $n.\n"; else $str .= "The matrix has in row $k unfortunately $len columns instead $n.\n"; return false; } } if ($m < 2 || $n < 2) { $str .= "At least two rows and at least two columns are required.\n"; return false; } for ($i = 0; $i < $m; $i++) for ($j = 0; $j < $n; $j++) $a[$i][$j] = floatval($a[$i][$j]); if (zero($a, $i, $j, $m, $n, $solve)) { return true; } else { for ($l = 1; $l < $n; $l++) $nbv[$l] = $l; for ($k = 1; $k < $m; $k++) $bv[$k] = $n + $k - 1; $str .= "The linear programme has $m rows and $n columns.\n"; trans($a, $i, $j, $m, $n, $bv, $nbv, $opt, $prog, $solve, $str); return false; } } else { $str .= "The input is either not sufficient or not feasible this way, because the numbers are e.g. not formed correctly.\n"; return false; } } function pivot(&$a, $cq, $cs, $d1, $d2, $d3, &$k, &$l, $n, $inc, &$max, &$min, $mem, $v, &$bet, $rule) { $ele = $a[0][$v]; switch ($rule) { case 1 : /* Compute pivot by steepest edge */ $quo = $ele * $ele * $d1[$v]; if ($quo > $max) { $k = $mem; $l = $v; $max = $quo; } elseif ($quo == $max) { $pro = $ele * $inc; if ($pro > $bet) { $k = $mem; $l = $v; $bet = $pro; } } break; case 2 : /* Compute pivot by best value */ $pro = $ele * $inc; if ($pro > $max) { $k = $mem; $l = $v; $max = $pro; } elseif ($pro == $max) { $quo = $ele * $ele * $d1[$v]; if ($quo > $bet) { $k = $mem; $l = $v; $bet = $quo; } } break; case 3 : /* Compute pivot by smallest angle */ $pro = 1 / $a[$mem][$v]; $per = $ele * $pro; $sum = $cs + $d1[$mem] * $pro - $per; $sca = $cq + ($d2[$mem] + $d3[$mem] * $pro) * $pro + $per * $per; $quo = abs($sum) * $sum / $sca; if ($quo < $min) { $k = $mem; $l = $v; $min = $quo; } elseif ($quo == $min) { $pro = $ele * $inc; if ($pro > $max) { $k = $mem; $l = $v; $max = $pro; } } break; } } function transform(&$a, $bv, &$mbv, &$cb, &$cn, &$i, &$j, $k, $l, $m, $n, $red, &$solve, $s, $val) { /* Compute right side */ $pivot = 1 / $a[$k][$l]; $a[$k][$l] = $pivot; $pline = $a[$k]; $rowel = -$pivot * $pline[0]; if ($val == 0) { for ($i = 0; $i < $k; $i++) $a[$i][0] += $rowel * $a[$i][$l]; for ($i = $k + 1; $i < $m; $i++) $a[$i][0] += $rowel * $a[$i][$l]; } else { for ($i = 0; $i < $val; $i++) { $r = $s[$i]; $a[$r][0] += $rowel * $a[$r][$l]; } } /* Compute pivot row */ $a[$k][0] = -$rowel; for ($j = 1; $j < $l; $j++) $a[$k][$j] *= $pivot; for ($j = $l + 1; $j < $n; $j++) $a[$k][$j] *= $pivot; $pivot = -$pivot; /* Compute pivot column and apply rectangle rule */ for ($i = 0; $i < $k; $i++) { $a[$i][$l] *= $pivot; $rowel = $a[$i][$l]; for ($j = 1; $j < $l; $j++) $a[$i][$j] += $rowel * $pline[$j]; for ($j = $l + 1; $j < $n; $j++) $a[$i][$j] += $rowel * $pline[$j]; } for ($i = $k + 1; $i < $m; $i++) { $a[$i][$l] *= $pivot; $rowel = $a[$i][$l]; for ($j = 1; $j < $l; $j++) $a[$i][$j] += $rowel * $pline[$j]; for ($j = $l + 1; $j < $n; $j++) $a[$i][$j] += $rowel * $pline[$j]; } zero($a, $i, $j, $m, $n, $solve); if ($red) redundant($a, $bv, $mbv, $cb, $cn, $m, $n); } function tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, $piv, $prog, $string, &$str) { /* Format tableau */ if ($opt == 1) $a[0][0] = -$a[0][0]; $str .= $string . '' . "\n"; $str .= ""; for ($j = 1; $j < $n; $j++) $str .= ""; $str .= "\n"; for ($i = 1; $i < $k; $i++) { $str .= ""; for ($j = 0; $j < $n; $j++) $str .= ""; $str .= "\n"; } if ($k > 0) { $str .= ""; for ($j = 0; $j < $l; $j++) $str .= ""; $str .= ($piv) ? '" : ""; for ($j = $l + 1; $j < $n; $j++) $str .= ""; $str .= "\n"; } for ($i = $k + 1; $i < $m; $i++) { $str .= ""; for ($j = 0; $j < $n; $j++) $str .= ""; $str .= "\n"; } if ($prog == 1) $str .= ($opt == 1 ) ? "" : ""; else $str .= ($opt == 1 ) ? "" : ""; for ($j = 0; $j < $l; $j++) $str .= ""; $str .= ($piv && $k == 0) ? '" : ""; for ($j = $l + 1; $j < $n; $j++) $str .= ""; $str .= "\n
BV:RS\\NBV:" . $nbv[$j] . "
" . $bv[$i] . "" . format($a[$i][$j]) . "
" . $bv[$k] . "" . format($a[$k][$j]) . "»' . format($a[$k][$l]) . "«" . format($a[$k][$l]) . "" . format($a[$k][$j]) . "
" . $bv[$i] . "" . format($a[$i][$j]) . "
PP Max:
PP Min:
DP Min:
DP Max:" . format($a[0][$j]) . "»' . format($a[0][$l]) . "«" . format($a[0][$l]) . "" . format($a[0][$j]) . "
\n"; if ($opt == 1) $a[0][0] = -$a[0][0]; } function loop(&$a, $m, $n, &$bv, &$nbv, &$mbv, &$cb, &$cn, $output, $red, $opt, &$perts, $prog, &$steps, &$str, &$solve, $rule, $last) { $cq = $cs = $l = $nul = 0; $d1 = array(); $d2 = array(); $d3 = array(); $d1[0] = $d2[0] = $d3[0] = 0; do { $k = $q = 0; /* Determine the coefficients of the objective function that are greater than zero and their indices */ for ($j = 1; $j < $n; $j++) if ($a[0][$j] > 0) $u[$q++] = $j; switch ($rule) { case 3 : /* Preparations smallest angle */ for ($i = 1; $i < $m; $i++) { $cq = 0; $cs = 0; $pro = 0; $quo = 0; $sca = 0; for ($j = 1; $j < $n; $j++) { $aij = $a[$i][$j]; $ele = $a[0][$j]; $mem = $ele * $ele; $inc = $mem * $aij; $cq += $mem; $cs += $ele; $pro += $inc; $quo += $inc * $aij; $sca += $ele * $aij; } $d1[$i] = -$sca; $d2[$i] = -$pro * 2; $d3[$i] = $quo; } break; default : /* Preparations steepest edge */ for ($j = 1; $j < $n; $j++) { $sca = 1; for ($i = 1; $i < $m; $i++) $sca += $a[$i][$j] * $a[$i][$j]; $d1[$j] = 1 / $sca; } } if ($q > 0) { $bet = $max = $o = $pro = 0; $min = INFINITY; /* Separate right sides that are zero from the other ones */ if ($nul == 0) { $p = 0; $s = null; for ($i = 1; $i < $m; $i++) if ($a[$i][0] == 0) $s[$o++] = $i; else $t[$p++] = $i; } else { for ($i = 1; $i < $m; $i++) if ($a[$i][0] == 0) $o++; } if ($o > 0 || $nul > 0) { if ($nul == 0) $nul = $o; /* Compute scalar products if need be */ if ($o > 0) { for ($i = 0; $i < $nul && $solve; $i++) { $sca = 0; $r = $s[$i]; for ($j = 1; $j < $n; $j++) $sca += $a[$r][$j] * $a[$r][$j]; $a[$r][0] += Sqrt($sca); if (!is_finite(floatval($a[$r][0]))) { trans($a, $r, 0, $m, $n, $bv, $nbv, $opt, $prog, $solve, $str); return $l; } } } /* Determine smallest step width if some right sides are zero */ for ($j = 0; $j < $q; $j++) { $v = $u[$j]; for ($i = 0; $i < $nul && $a[$s[$i]][$v] <= 0; $i++); if ($i == $nul) { $inc = INFINITY; for ($i = 0; $i < $p; $i++) { $r = $t[$i]; $ele = $a[$r][$v]; if ($ele > 0) { $quo = $a[$r][0] / $ele; if ($quo < $inc) { $inc = $quo; $mem = $r; } } } if ($inc < INFINITY) pivot($a, $cq, $cs, $d1, $d2, $d3, $k, $l, $n, $inc, $max, $min, $mem, $v, $bet, $rule); } } /* Undo perturbations */ if ($k > 0) { for ($i = 0; $i < $nul; $i++) $a[$s[$i]][0] = 0; $nul = 0; } } /* Determine smallest step width if all right sides are not zero */ else { for ($j = 0; $j < $q; $j++) { $v = $u[$j]; $inc = INFINITY; for ($i = 1; $i < $m; $i++) { $ele = $a[$i][$v]; if ($ele > 0) { $quo = $a[$i][0] / $ele; if ($quo < $inc) { $inc = $quo; $mem = $i; } } } if ($inc < INFINITY) pivot($a, $cq, $cs, $d1, $d2, $d3, $k, $l, $n, $inc, $max, $min, $mem, $v, $bet, $rule); } } /* Determine smallest step width if the objective function could not be increased */ if ($nul > 0) { $bet = $max = 0; $min = INFINITY; for ($j = 0; $j < $q; $j++) { $v = $u[$j]; $inc = INFINITY; for ($i = 0; $i < $nul; $i++) { $r = $s[$i]; $ele = $a[$r][$v]; if ($ele > 0) { $quo = $a[$r][0] / $ele; if ($quo < $inc) { $inc = $quo; $mem = $r; } } } if ($inc < INFINITY) pivot($a, $cq, $cs, $d1, $d2, $d3, $k, $l, $n, $inc, $max, $min, $mem, $v, $bet, $rule); } $per = ($k > 0) ? ++$perts : false; } else { $per = false; } /* Transform tableau and build details if demanded */ if ($k > 0) { swap($bv[$k], $nbv[$l]); $steps++; if ($output) { tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); $str .= "
Das Pivotelement ist a[$k][" . ($l + 1) . "] = " . format($a[$k][$l]) . ".
\n"; if ($per) $str .= "
Eine Perturbation war erforderlich.
\n"; } $mem = $a[0][0]; transform($a, $bv, $mbv, $cb, $cn, $i, $j, $k, $l, $m, $n, $red, $solve, $s, $nul); if ($last && $a[0][0] != $mem) return; if (!$solve) { trans($a, $i, $j, $m, $n, $bv, $nbv, $opt, $prog, $solve, $str); $k = 0; } } } } while ($k > 0); /* Undo perturbations */ if ($nul > 0 && $solve) { if ($output) tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, false, $prog, "", $str); for ($i = 0; $i < $nul; $i++) $a[$s[$i]][0] = 0; if ($output) $str .= "
Tableau without perturbations:
\n"; } return $l; } function numbers($perts, $steps, $string, &$str) { /* Prepare pivot steps and perturbations for output */ if ($steps == 1) $str .= "
There were$string one pivot step and "; else $str .= "
There were$string $steps pivot steps and "; if ($perts == 0) $str .= "no perturbations needed.
\n"; else $str .= ($perts == 1) ? "one perturbation needed.
\n" : "$perts perturbations needed.\n"; } function results(&$a, &$e, &$n1, &$n2, $m, $n, &$bv, &$nbv, &$mbv, $cb, $cn, $red, $opt, $perts, $prog, $steps, &$str, &$solve, $nor1, $scal, $nor2, $rule) { /* Collect results */ numbers($perts, $steps, '', $str); if ($prog == 1) $str .= ($opt == 1) ? "
The maximum is " . format(-$a[0][0]) . "." : "
The minimum is " . format($a[0][0]) . "."; else $str .= ($opt == 1) ? "
The minimum is " . format(-$a[0][0]) . "." : "
The maximum is " . format($a[0][0]) . "."; $x = array_fill(1, $n - 1, 0); $y = array_fill(1, $m - 1, 0); for ($i = 1; $i < $m; $i++) if ($bv[$i] < $n) $x[$bv[$i]] = $a[$i][0]; if ($scal) for ($j = 1; $j < $n; $j++) $x[$j] *= $e[$j]; $str .= ($prog == 1) ? "
\nPrimal solution: " : "
\nDual   solution: "; for ($j = 1; $j < $n; $j++) { $str .= $x[$j] . " "; if ($nbv[$j] >= $n) $y[$nbv[$j] - $n + 1] = -$a[0][$j]; } if ($nor1 || $nor2) for ($i = 1; $i < $m; $i++) $y[$i] *= $n1[$i]; if ($nor1 && $scal && $nor2) for ($i = 1; $i < $m; $i++) $y[$i] *= $n2[$i]; $str .= ($prog == 1) ? "
\nDual   solution: " : "
\nPrimal solution: "; for($i = 1; $i < $m; $i++) $str .= $y[$i] . " "; $str .= "
\n"; other($a, $m, $n, $bv, $nbv, $mbv, $cb, $cn, $red, $opt, $perts, $prog, $steps, $str, $solve, $rule); } function other(&$a, $m, $n, &$bv, &$nbv, &$mbv, $cb, $cn, $red, $opt, &$perts, $prog, &$steps, &$str, &$solve, $rule) { /* Check if there are other solutions */ $q = 0; for ($j = 1; $j < $n; $j++) if ($a[0][$j] == 0) $u[$q++] = $j; if ($q > 0) { $none = true; for ($j = 0; $j < $q; $j++) { $v = $u[$j]; for ($i = 1; $i < $m && ($a[$i][$v] <= 0 && $a[$i][0] == 0 || $a[$i][0] != 0); $i++); $w[$j] = $i == $m; } for ($j = 0; $j < $q && $none; $j++) if ($nbv[$u[$j]] < $n && $w[$j]) $none = false; for ($i = 1; $i < $m && $none; $i++) if ($bv[$i] < $n) for ($j = 0; $j < $q; $j++) if ($a[$i][$u[$j]] != 0 && $w[$j]) $none = false; if ($none) { if ($q > 1) { $l = 0; for ($j = 1; $j < $n; $j++) { if ($a[0][$j] == 0) { if ($j > ++$l) { for ($i = 1; $i < $m; $i++) $a[$i][$l] = $a[$i][$j]; $a[0][$l] = 1; } else { $a[0][$j] = 1; } } } $n = ++$l; for ($j = 1; $j < $n; $j++) $nbv[$j] = $j; for ($i = 1; $i < $m; $i++) $bv[$i] = $n + $i - 1; $mem = $a[0][0]; loop($a, $m, $n, $bv, $nbv, $mbv, $cb, $cn, false, false, $opt, $perts, $prog, $steps, $str, $solve, $rule, true); if ($solve) { $i = 0; if ($a[0][0] == $mem) { $q = 0; for ($j = 1; $j < $n; $j++) if ($a[0][$j] >= 0) $u[$q++] = $j; if ($q > 0) { for ($j = 0; $j < $q && $i < $m; $j++) { $v = $u[$j]; for ($i = 1; $i < $m && $a[$i][$v] <= 0; $i++); } } } $str .= ($a[0][0] != $mem || $i == $m) ? (($prog == 1) ? "
There are further primal solutions.
\n" : "
There are further dual solutions.
\n") : (($prog == 1) ? "
There are no further primal solutions.
\n" : "
There are no further dual solutions.
\n"); } } else { $str .= ($prog == 1) ? "
There are no further primal solutions.
\n" : "
There are no further dual solutions.
\n"; } } else { $str .= ($prog == 1) ? "
There are further primal solutions.
\n" : "
There are further dual solutions.
\n"; } } else { $str .= ($prog == 1) ? "
There are no further primal solutions.
\n" : "
There are no further dual solutions.
\n"; } } /* Read file, prepare it for output with number of rows and columns, build main array */ $fp = fopen($file, 'r+'); $str .= "
Content of the file $file:\n
"; $lines = ''; while ($line = fgets($fp)) { if ($line) $lines .= "$line
"; } fclose($fp); $a = array(); $bv = array(); $nbv = array(); $n1 = array(); $n2 = array(); $m = $n = 0; $solve = true; $str .= "$lines\n"; if (check($a, $m, $n, $bv, $nbv, $opt, $prog, $solve, $str, $lines)) { /* Build tableau, dualise if requested and index (non-) base variables */ if ($opt == 1) $a[0][0] = -$a[0][0]; else for ($j = 1; $j < $n; $j++) $a[0][$j] = -$a[0][$j]; if ($op == 2) for ($i = 1; $i < $m; $i++) for ($j = 0; $j < $n; $j++) $a[$i][$j] = -$a[$i][$j]; if ($prog == 2) { $opt = ($opt == 1) ? 2 : 1; $min = min($m, $n); for ($i = 0; $i < $min; $i++) { for ($j = 0; $j <= $i; $j++) { $tmp = $a[$i][$j]; $a[$i][$j] = -$a[$j][$i]; $a[$j][$i] = -$tmp; } } if ($m != $n) { if ($m < $n) { for ($j = $m; $j < $n; $j++) { for ($i = 0; $i < $m; $i++) { $a[$j][$i] = -$a[$i][$j]; } } } else { for ($i = $n; $i < $m; $i++) { for ($j = 0; $j < $n; $j++) { $a[$j][$i] = -$a[$i][$j]; } } } swap($m, $n); } } for ($j = 1; $j < $n; $j++) $nbv[$j] = $j; for ($i = 1; $i < $m; $i++) $bv[$i] = $n + $i - 1; $str .= "The linear programme has $m rows and $n columns.
\n"; /* Check whether the linear programme can be solved for simple constraints */ for ($i = 1; $i < $m && $solve; $i++) { for ($j = 1; $j < $n && $a[$i][$j] >= 0; $j++); if ($j == $n) { if ($a[$i][0] < 0) { $solve = false; tableau($a, $i, 0, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); $str .= "
The linear programme cannot be solved.
\n"; } else { if ($a[$i][0] == 0) { for ($j = 1; $j < $n && $a[$i][$j] > 0; $j++); if ($j == $n) { $solve = false; tableau($a, $i, 0, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); for ($i = 1; $i < $m && $a[$i][0] >= 0; $i++); if ($i == $m) { if ($prog == 1) $str .= ($opt == 1) ? "
The maximum is " . format(-$a[0][0]) . "." : "
The minimum is " . format($a[0][0]) . "."; else $str .= ($opt == 1) ? "
The minimum is " . format(-$a[0][0]) . "." : "
The maximum is " . format($a[0][0]) . "."; $str .= ($prog == 1) ? "
\nPrimal solution: " : "
\nDual solution: "; for ($j = 1; $j < $n; $j++) $str .= "0 "; $str .= ($prog == 1) ? "
\nDual solution: " : "
\nPrimal solution: "; for ($i = 1; $i < $m; $i++) $str .= "0 "; $str .= "
\n"; } else { $str .= "
The linear programme cannot be solved.
\n"; } } } } } } /* Initialise variables, put normalised rows into new matrix and determine trivial redundant constraints if any */ $phase1 = false; if ($red) { $cb = 0; $cn = 0; $mbv = array(); $mbv = array_fill(0, $m + $n - 1, true); $b = $a; for ($i = 1; $i < $m; $i++) { $max = $sca = 0; for ($j = 1; $j < $n; $j++) { $abs = abs($b[$i][$j]); if ($abs > $max) $max = $abs; } $max = ($max == 0) ? 1 : 1 / $max; $b[$i][0] *= $max; for ($j = 1; $j < $n; $j++) { $b[$i][$j] *= $max; $sca += $b[$i][$j] * $b[$i][$j]; } $sca = ($sca > 0) ? 1 / sqrt($sca) : 1; for ($j = 0; $j < $n; $j++) $b[$i][$j] *= $sca; $f[$i] = $b[$i][0] >= 0; if ($f[$i]) { for ($j = 1; $j < $n && $b[$i][$j] <= 0; $j++); $f[$i] = $j == $n; } } /* Compare rows and compute upper and lower limits to remove further redundant rows if possible */ if ($m > 2) { for ($i = 1; $i < $m; $i++) { $hsi = $b[$i][0]; $hgi = $hsi > 0; $hli = $hsi < 0; $hge = $hsi >= 0; $hle = $hsi <= 0; $heq = $hsi == 0; for ($k = $i + 1; !$f[$i] && $k < $m; $k++) { if (!$f[$k]) { $lsi = $b[$k][0]; if ($hle && $lsi >= 0) $l = $i; else if ($hge && $lsi <= 0) $l = $k; else $l = 0; $h = ($l == $i) ? $k : $i; if ($l != 0) { $llim = 0; $ulim = INFINITY; } else { $l = $k; if ($hgi) { $llim = 0; $ulim = $hsi / $lsi; } else { $llim = $hsi / $lsi; $ulim = INFINITY; } } $one = true; for ($j = 1; $one && $j < $n; $j++) { $hco = $b[$h][$j]; $lco = $b[$l][$j]; $hgo = $hco > 0; $hlo = $hco < 0; $lgo = $lco > 0; $llo = $lco < 0; if ($hgo && $lgo || $hlo && $llo) { $hco /= $lco; if ($lgo && $llim < $hco) $llim = $hco; elseif ($llo && $ulim > $hco) $ulim = $hco; if ($llim > $ulim + EPSILON) $one = false; } else { if ($hgo && $lco <= 0) $one = false; elseif ($hco >= 0 && $llo) $one = false; } } $f[$h] = $one; /* Alter probable redundant row if necessary and compute again */ if (!$one && ($hgi && $lsi > 0 || $hli && $lsi < 0 || $heq && $lsi == 0)) { $llim = ($hli) ? $lsi / $hsi : 0; $ulim = ($hgi) ? $lsi / $hsi : INFINITY; $one = true; for ($j = 1; $one && $j < $n; $j++) { $hco = $b[$h][$j]; $lco = $b[$l][$j]; $hgo = $hco > 0; $hlo = $hco < 0; $lgo = $lco > 0; $llo = $lco < 0; if ($hgo && $lgo || $hlo && $llo) { $lco /= $hco; if ($hgo && $llim < $lco) $llim = $lco; elseif ($hlo && $ulim > $lco) $ulim = $lco; if ($llim > $ulim + EPSILON) $one = false; } else { if ($hlo && $lco >= 0) $one = false; elseif ($hco <= 0 && $lgo) $one = false; } } $f[$l] = $one; } } } } } /* Output redundant rows if any and remove them from the initial matrix */ $q = 0; for ($i = 1; $i < $m; $i++) if ($f[$i]) $g[$q++] = $i; if ($q > 0) { $str .= ($q == 1) ? "
The following row is" : "
The following rows are"; $str .= " redundant:
\n" . '' . "\n"; for ($i = 0; $i < $q; $i++) { $str .= ""; for ($j = 0; $j < $n; $j++) $str .= ""; $str .= "\n"; } $str .= "
" . $a[$g[$i]][$j] . "
\n"; $q = 1; for ($i = 1; $i < $m; $i++) { if (!$f[$i]) { for ($j = 0; $j < $n; $j++) $a[$q][$j] = $a[$i][$j]; $q++; } } $m = $q; } $b = $a; } if ($solve) { /* Check whether linear programme is limited */ $limit = true; for ($j = 1; $j < $n && $limit; $j++) { if ($a[0][$j] > 0) { for ($i = 1; $i < $m && $a[$i][$j] < 0; $i++); $limit = $i < $m; } } if ($limit && $m > 1) { /* Initialise variables and determine smallest negative value of the right side if any */ $min = $perts = $steps = 0; for ($i = 1; $i < $m; $i++) { $ele = $a[$i][0]; if ($ele < $min) { $k = $i; $min = $ele; } } if ($min < 0) { /* Phase I: Add slack variable in the last column, save objective function and make tableau valid by pivot step */ $phase1 = true; if ($output) tableau($a, $k, 0, $m, $n, $bv, $nbv, $opt, false, $prog, "", $str); /* Normalise and scale if demanded */ if ($nor1 || $nor2) { if ($output) $str .= "
The 1st normalisation will be executed.
\n"; normalise($a, $n1, $i, $j, $m, $n); } if ($scal) { if ($output) $str .= "
The scaling will be executed.
\n"; scale($a, $e, $i, $j, $m, $n); if ($nor1 && $nor2) { if ($output) $str .= "
The 2nd normalisation will be executed.
\n"; normalise($a, $n2, $i, $j, $m, $n); } } for ($i = 1; $i < $m; $i++) $a[$i][$n] = -1; $c = $a[0]; $mem = $c[0]; $a[0] = array_fill(0, $n, 0); $a[0][$n] = -1; $nbv[$n] = $bv[$k]; $bv[$k] = $c[0] = 0; $none = true; for ($j = 1; $j < $n; $j++) $none &= $c[$j] == 0; $l = $n++; $steps++; if ($output) { tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, true, $prog, "
Phase I:
\n", $str); $str .= "
The pivot is a[$k][" . ($l + 1) . "] = -1.
\n"; } transform($a, $bv, $mbv, $cb, $cn, $i, $j, $k, $l, $m, $n, $red, $solve, null, 0); if ($solve) { /* Loop, check feasibility and collect results */ $l = loop($a, $m, $n, $bv, $nbv, $mbv, $cb, $cn, $output, $red, $opt, $perts, $prog, $steps, $str, $solve, $rule, false); if ($solve) { if ($a[0][0] != 0) { $solve = false; // infeasible tableau($a, 0, 0, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); $str .= "
The linear programme cannot be solved.
\n"; numbers($perts, $steps, '', $str); } else { tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, false, $prog, "", $str); if ($nbv[$l] != 0) { /* Move slack variable from the base variables to the last column of the non-base variables and drop last column */ for ($l = $n - 1; $l > 0 && $a[$k][$l] == 0; $l--); $bv[$k] = $nbv[$l]; $nbv[$l] = 0; $steps++; transform($a, $bv, $mbv, $cb, $cn, $i, $j, $k, $l, $m, $n, $red, $solve, null, 0); if ($output) tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, true, $prog, "
The pivot is a[$k][" . ($l + 1) . "] = " . format(1 / $a[$k][$l]) . ".
\n", $str); } if ($solve) { $phase1 = false; $n--; if ($nbv[$l] == 0 && $l < $n) { /* Replace column in tableau with slack variable by last column and drop the latter */ for ($i = 1; $i < $m; $i++) $a[$i][$l] = $a[$i][$n]; $nbv[$l] = $nbv[$n]; } $a[0][0] = $mem; if ($none) { /* All coefficients of the objective function are zero */ $l = 0; tableau($a, $k, $l, $m, $n, $bv, $nbv, $opt, false, $prog, "
End of phase I:
\n", $str); results($a, $e, $n1, $n2, $m, $n, $bv, $nbv, $mbv, $cb, $cn, $red, $opt, $perts, $prog, $steps, $str, $solve, $nor1, $scal, $nor2, $rule); $solve = false; } else { /* Translate objective function into new variables */ for ($j = 1; $j < $n; $j++) { $l = $nbv[$j]; $a[0][$j] = ($l < $n) ? $c[$l] : 0; } for ($i = 1; $i < $m; $i++) { $l = $bv[$i]; if ($l < $n) { $obj = -$c[$l]; for ($j = 0; $j < $n; $j++) $a[0][$j] += $obj * $a[$i][$j]; } } for ($j = 0; $j < $n; $j++) if (abs($a[0][$j]) < EPSILON) $a[0][$j] = 0; if ($red) redundant($a, $bv, $mbv, $cb, $cn, $m, $n); if ($output) $str .= "
End of phase I:
\n"; } } else { $str .= "
The boundary of computation was transcended.
\n"; numbers($perts, $steps, '', $str); } } } else { numbers($perts, $steps, '', $str); } } else { $str .= "
The boundary of computation was transcended.
\n"; numbers($perts, $steps, '', $str); } } else { if ($output) tableau($a, 0, 0, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); /* Normalise and scale if demanded */ if ($nor1 || $nor2) { if ($output) $str .= "
The 1st normalisation will be executed.
\n"; normalise($a, $n1, $i, $j, $m, $n); } if ($scal) { if ($output) $str .= "
The scaling will be executed.
\n"; scale($a, $e, $i, $j, $m, $n); if ($nor1 && $nor2) { if ($output) $str .= "
The 2nd normalisation will be executed.
\n"; normalise($a, $n2, $i, $j, $m, $n); } } } /* Phase II: Loop, check whether the linear programme is limited and collect further results */ if ($solve) { $phase1 = false; for ($j = 1; $j < $n && $a[0][$j] <= 0; $j++); if ($j < $n) { if ($min < 0) { if ($output) { numbers($perts, $steps, '', $str); $str .= "
Phase II:
\n"; } $allp = $perts; $alls = $steps; $perts = $steps = 0; } loop($a, $m, $n, $bv, $nbv, $mbv, $cb, $cn, $output, $red, $opt, $perts, $prog, $steps, $str, $solve, $rule, false); if ($solve) { /* Check whether linear programme is limited */ for ($j = 1; $j < $n && $a[0][$j] <= 0; $j++); if ($j < $n) { tableau($a, 0, $j, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); $j++; $solve = false; $str .= "
The linear programme is unlimited in column $j.
\n"; numbers($perts, $steps, '', $str); } else { tableau($a, 0, 0, $m, $n, $bv, $nbv, $opt, false, $prog, "", $str); results($a, $e, $n1, $n2, $m, $n, $bv, $nbv, $mbv, $cb, $cn, $red, $opt, $perts, $prog, $steps, $str, $solve, $nor1, $scal, $nor2, $rule); } } else { numbers($perts, $steps, '', $str); } if ($min < 0) numbers($allp + $perts, $alls + $steps, ' in total', $str); } else { tableau($a, 0, 0, $m, $n, $bv, $nbv, $opt, false, $prog, "", $str); results($a, $e, $n1, $n2, $m, $n, $bv, $nbv, $mbv, $cb, $cn, $red, $opt, $perts, $prog, $steps, $str, $solve, $nor1, $scal, $nor2, $rule); } } } else { if ($limit) { tableau($a, 0, 0, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); if ($prog == 1) $str .= ($opt == 1) ? "
The maximum is " . format(-$a[0][0]) . "." : "
The minimum is " . format($a[0][0]) . "."; else $str .= ($opt == 1) ? "
The minimum is " . format(-$a[0][0]) . "." : "
The maximum is " . format($a[0][0]) . "."; $str .= ($prog == 1) ? "
\nPrimale Lösung: " : "
\nDuale Lösung: "; for ($j = 1; $j < $n; $j++) $str .= "0 "; $str .= ($prog == 1) ? "
\nDual solution: Does not exist.
\nThere are no further primal solutions.
\n" : "
\nPrimal solution: Does not exist.
\nThere are no further dual solutions.
\n"; } else { tableau($a, 0, $j - 1, $m, $n, $bv, $nbv, $opt, true, $prog, "", $str); $str .= "
The linear programme is unlimited in column $j.
\n"; } } } if ($red) { /* Output redundant rows resp. nonnegativity conditions */ if ($phase1) $n--; if ($cb > 0) { $too = ($q > 0) ? " also" : ""; $str .= ($cb == 1) ? "
The following row is" : "
The following rows are"; $str .= "$too redundant:
\n" . '' . "\n"; for ($i = 1; $i < $m; $i++) { if (!$mbv[$i + $n - 1]) { $str .= ""; for ($j = 0; $j < $n; $j++) $str .= ""; $str .= "\n"; } } $str .= "
" . $b[$i][$j] . "
\n"; } if (!$mbv[0]) { $cn--; $mbv[0] = true; } if ($cn > 0) { $str .= ($cn == 1) ? "
The following variable doesn't" : "
The following variables don't"; $str .= " need a nonnegativity condition:
\n"; for ($j = 1; $j < $n; $j++) if (!$mbv[$j]) $str .= $j + " "; $str .= "
\n"; } } } $str .= "

© 2018 by Boris Haase
\n"; $str .= 'top

' . "\n\n"; /* Output */ echo $str; ?>