Симплекс-метод с равенствами

    /// <summary>
    /// Система линейных уравнений/неравенств
    /// </summary>
    public class LinearEquationSystem : List<LinearEquation>
    {
        public HashSet<string> Variables { get; } = new HashSet<string>();
 
        public new void Add(LinearEquation eq)
        {
            foreach (var v in eq.Keys)
                Variables.Add(v);
 
            base.Add(eq);
        }
 
        //парсим систему уравнений
        public static LinearEquationSystem Parse(string str)
        {
            var res = new LinearEquationSystem();
 
            foreach (var line in str.Split('\r', '\n').Where(s => !string.IsNullOrWhiteSpace(s)))
                res.Add(LinearEquation.Parse(line));
 
            return res;
        }
    }

    /// <summary>
    /// Линейное уравнение/неравенство.
    /// Содержит пары переменная-коэффициент.
    /// Свободный члнен хранится как перменная с пустым именем.
    /// Уравнение автоматически приводится к каноническому виду - слева перменные и свободный член, справа - ноль.
    /// Одноименные перменные приводятся к одному члену.
    /// </summary>
    [DebuggerDisplay("{ToString()}")]
    public class LinearEquation : Dictionary<string, double>
    {
        public string Sign { get; private set; }
 
        //добавление члена уравнения
        public new void Add(string varName, double coeff)
        {
            this[varName] = GetCoeff(varName) + coeff;
        }
 
        public double GetCoeff(string varName)
        {
            var c = 0d;
            TryGetValue(varName, out c);
            return c;
        }
 
        //парсим уравнение
        public static LinearEquation Parse(string str)
        {
            var result = new LinearEquation();
            var m = Regex.Match(str, @"^([^\!<>=]+?)([\!<>=]+)([^\!<>=]+?)$");
            if (!m.Success)
                throw new Exception("Syntax error");
            Parse(result, m.Groups[1].Value, 1);//left part
            Parse(result, m.Groups[3].Value, -1);//right part
            result.Sign = m.Groups[2].Value;
 
            return result;
        }
 
        //парсим выражение
        private static void Parse(LinearEquation result, string str, int coeff)
        {
            foreach (Match m in Regex.Matches(str, @"([+\-]?)\s*([\w*\.\s]+)"))
            {
                var s = m.Groups[2].Value.Trim();
 
                if (s == "") continue;
 
                var v = s;
                var c = 1d;
 
                //выделяем коэффициент и имя переменной
                var cm = Regex.Match(s, @"(?<coeff>\d+)\s*\*\s*(?<var>[a-z0-9]*)|(?<coeff>\d+)(?<var>[a-z0-9]*)|(?<var>[a-z0-9][a-z0-9]*)");
                if (cm.Success && cm.Groups["var"].Value != "")
                {
                    v = cm.Groups["var"].Value;
                    if(cm.Groups["coeff"].Value != "")
                        c = double.Parse(cm.Groups["coeff"].Value, CultureInfo.InvariantCulture);
                }
                else
                {
                    if (double.TryParse(s, NumberStyles.AllowDecimalPoint, CultureInfo.InvariantCulture, out c))
                        v = "";
                    else
                        c = 1d;
                }
 
                //обрабатываем знак
                if (m.Groups[1].Value == "-")
                    c *= -1;
 
                c *= coeff;
 
                result.Add(v, c);
            }
        }
 
        public override string ToString()
        {
            var sb = new StringBuilder();
            foreach (var pair in this)
            {
                var mult = string.IsNullOrEmpty(pair.Key) ? "" : "*";
                sb.AppendFormat("{0}{1:G}{3}{2}", Math.Sign(pair.Value) == -1 ? "-" : "+", Math.Abs(pair.Value), pair.Key, mult);
            }
 
            sb.AppendFormat(" {0} 0", Sign);
 
            return sb.ToString().Trim('+');
        }
 
        public void Invert()
        {
            switch (Sign)
            {
                case ">=": Sign = "<="; break;
                case "<=": Sign = ">="; break;
                case ">": Sign = "<"; break;
                case "<": Sign = ">"; break;
                case "=>": Sign = "=<"; break;
            }
 
            foreach (var v in Keys.ToArray())
                this[v] *= -1;
        }
    }

    public class SimplexMethod
    {
        public const double Eps = 1e-9;
 
        public SimplexMethodResult Solve(double[,] A, double[] b, double[] c)
        {
            var n = A.GetLength(1);
            var m = A.GetLength(0);
 
            var ANew = new double[m,n + m];
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    ANew[i, j] = A[i, j];
                }
            }
            A = ANew;
            var cNew = new double[n + m];
            for (int i = 0; i < n; i++)
            {
                cNew[i] = c[i];
            }
            c = cNew;
            n = n + m;
 
            var N = new HashSet<int>();
            for (int i = 0; i < n - m; i++)
            {
                N.Add(i);
            }
            var B = new HashSet<int>();
            var id = new int[n + 1];
            for (int i = n - m; i < n; i++)
            {
                B.Add(i);
                id[i] = i - (n - m);
            }
 
            int l;
            int e;
            double v = 0;
            if (!b.All(value => value > -Eps))
            {
                ANew = new double[m, n + 1];
                for (int i = 0; i < m; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        ANew[i, j] = A[i, j];
                    }
                    ANew[i, n] = -1;
                }
                A = ANew;
 
                cNew = new double[n + 1];
                cNew[n] = -1;
 
                N.Add(n);
 
                e = n;
                l = -1;
                foreach (var i in B)
                {
                    if (b[id[i]] < -Eps && (l == -1 || b[id[l]] > b[id[i]]))
                    {
                        l = i;
                    }
                }
 
                v = 0;
                var result = Pivot(N, B, A, id, b, cNew, ref v, l, e);
                A = result.A;
                b = result.b;
                cNew = result.c;
 
                while (N.Any(i => cNew[i] > Eps))
                {
                    e = -1;
                    foreach (var i in N)
                    {
                        if (cNew[i] > Eps)
                        {
                            e = i;
                            break;
                        }
                    }
 
                    l = -1;
                    foreach (var i in B)
                    {
                        if (A[id[i], e] > Eps && (l == -1 || b[id[l]] / A[id[l], e] > b[id[i]] / A[id[i], e]))
                        {
                            l = i;
                        }
                    }
 
                    if (l == -1)
                    {
                        return new SimplexMethodResult { Status = SimplexMethodResultStatus.Unlimited };
                    }
                    result = Pivot(N, B, A, id, b, cNew, ref v, l, e);
                    A = result.A;
                    b = result.b;
                    cNew = result.c;
                }
 
                if (v < - Eps)
                {
                    return new SimplexMethodResult { Status = SimplexMethodResultStatus.NoSolution };
                }
 
                if (B.Contains(n))
                {
                    B.Remove(n);
                    e = -1;
                    foreach (var j in N)
                    {
                        if (Math.Abs(A[id[n], j]) > Eps)
                        {
                            e = j;
                            break;
                        }
                    }
                    var d = A[id[n], e];
                    for (int j = 0; j < n; j++)
                    {
                        A[id[n], j] /= d;
                    }
                    A[id[n], e] = 0;
                    foreach (var i in B)
                    {
                        d = -A[id[i], e];
                        for (int j = 0; j < n; j++)
                        {
                            A[id[i], j] += d * A[id[n], j];
                        }
                        A[id[i], e] = 0;
                    }
                    N.Remove(e);
                    B.Add(e);
                    id[e] = id[n];
                }
                else if (N.Contains(n))
                {
                    N.Remove(n);
                }
 
                cNew = new double[n];
                for (int i = 0; i < n; i++)
                {
                    if (!B.Contains(i))
                    {
                        cNew[i] += c[i];
                    }
                    else
                    {
                        for (int j = 0; j < n; j++)
                        {
                            cNew[j] -= c[i] * A[id[i], j];
                        }
                        v += c[i] * b[id[i]];
                    }
                }
                c = cNew;
 
                ANew = new double[m, n];
                for (int i = 0; i < m; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        ANew[i, j] = A[i, j];
                    }
                }
                A = ANew;
            }
 
            while (N.Any(i => c[i] > Eps))
            {
                e = -1;
                foreach (var i in N)
                {
                    if (c[i] > Eps)
                    {
                        e = i;
                        break;
                    }
                }
 
                l = -1;
                foreach (var i in B)
                {
                    if (A[id[i], e] > Eps && (l == -1 || b[id[l]] / A[id[l], e] > b[id[i]] / A[id[i], e]))
                    {
                        l = i;
                    }
                }
 
                if (l == -1)
                {
                    return new SimplexMethodResult { Status = SimplexMethodResultStatus.Unlimited };
                }
                var result = Pivot(N, B, A, id, b, c, ref v, l, e);
                A = result.A;
                b = result.b;
                c = result.c;
            }
 
            var x = new double[n - m];
            for (int i = 0; i < n - m; i++)
            {
                if (B.Contains(i))
                {
                    x[i] = b[id[i]];
                }
            }
 
            return new SimplexMethodResult { Status = SimplexMethodResultStatus.Ok, Value = v, Solution = x };
        }
 
        private PivotResult Pivot(ISet<int> N, ISet<int> B, double[,] A, IList<int> id, IList<double> b, IList<double> c, ref double v, int l, int e)
        {
            var ANew = new double[A.GetLength(0), A.GetLength(1)];
            var bNew = new double[b.Count];
            var cNew = new double[c.Count];
 
            id[e] = id[l];
 
            bNew[id[e]] = b[id[l]] / A[id[l], e];
 
            foreach (var j in N)
            {
                if (j != e)
                {
                    ANew[id[e], j] = A[id[l], j] / A[id[l], e];
                }
            }
            ANew[id[e], l] = 1 / A[id[l], e];
 
            foreach (var i in B)
            {
                if (i != l)
                {
                    bNew[id[i]] = b[id[i]] - A[id[i], e] * bNew[id[e]];
                    foreach (var j in N)
                    {
                        if (j != e)
                        {
                            ANew[id[i], j] = A[id[i], j] - A[id[i], e] * ANew[id[e], j];
                        }
                    }
                    ANew[id[i], l] = -A[id[i], e] * ANew[id[e], l];
                }
            }
 
            v += c[e] * bNew[id[e]];
            foreach (var j in N)
            {
                if (j != e)
                {
                    cNew[j] = c[j] - c[e] * ANew[id[e], j];
                }
            }
            cNew[l] = -c[e] * ANew[id[e], l];
 
            N.Remove(e);
            N.Add(l);
            B.Remove(l);
            B.Add(e);
 
            return new PivotResult { A = ANew, c = cNew, b = bNew };
        }
    }
 
    public class SimplexMethodResult
    {
        public SimplexMethodResultStatus Status { get; set; }
 
        public double Value { get; set; }
 
        public double[] Solution { get; set; }
    }
 
    public enum SimplexMethodResultStatus
    {
        Ok,
        NoSolution,
        Unlimited
    }
 
    internal class PivotResult
    {
        public double[,] A { get; set; }
 
        public double[] c { get; set; }
 
        public double[] b { get; set; }
    }

    /// <summary>
    /// Парсит выражения, подготавливает матрицы для симплекс метода, решает симплекс
    /// </summary>
    public class SimplexSolver
    {
        public List<string> GetVarNames(LinearEquationSystem les)
        {
            var vars = les.Variables.Where(s => s != "max" && s != "min" && s != "").ToList();
            return vars;            
        }
 
        public void PrepareMatrix(IList<LinearEquation> les, List<string> vars, out double[,] A, out double[] b)
        {
            var cols = vars.Count;
            var rows = les.Count;
            A = new double[rows, cols];
            b = new double[rows];
            
            for (int j = 0; j < rows; j++)
            {
                var eq = les[j];
                b[j] = -eq.GetCoeff("");//свободный член
                for (int i = 0; i < cols; i++)
                    A[j, i] = les[j].GetCoeff(vars[i]);
            }
        }
 
        public void PrepareTarget(LinearEquation eq, List<string> vars, out double[] c)
        {
            var cols = vars.Count;
            c = new double[cols];
 
            for (int i = 0; i < cols; i++)
                c[i] = eq.GetCoeff(vars[i]);
        }
 
        public Dictionary<string, double> Solve(LinearEquationSystem les)
        {
            //получакем имена всех перменных
            var vars = GetVarNames(les);
 
            //подготавливаем матрицы
            double[,] A;//коэфф ограничений
            double[] b;//свободный член ограничений
            double[] c;//коэфф целевой функции
            
            PrepareMatrix(les.Where(e=>e.Sign == "=").ToList(), vars, out A, out b);
            PrepareTarget(les.FirstOrDefault(e => e.Sign == "=>" || e.Sign == "=<"), vars, out c);
 
            //решаем симплекс
            var simplex = new SimplexMethod();
            var solution = simplex.Solve(A, b, c);
 
            //подготавливаем результат
            if (solution.Status == SimplexMethodResultStatus.Ok)
            {
                var res = new Dictionary<string, double>();
                for (int i = 0; i < vars.Count; i++)
                    res[vars[i]] = solution.Solution[i];
 
                return res;
            }
 
            return null;
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            var task = @"
1*x000 + 1*x001 + 1*x002 + 1*x010 + 1*x011 + 1*x012 + 1*x020 + 1*x021 + 1*x022 = 1
1*x100 + 1*x101 + 1*x102 + 1*x110 + 1*x111 + 1*x112 + 1*x120 + 1*x121 + 1*x122 = 1
1*x200 + 1*x201 + 1*x202 + 1*x210 + 1*x211 + 1*x212 + 1*x220 + 1*x221 + 1*x222 = 1
1*x000 + 1*x001 + 1*x002 + 1*x100 + 1*x101 + 1*x102 + 1*x200 + 1*x201 + 1*x202 = 1
1*x010 + 1*x011 + 1*x012 + 1*x110 + 1*x111 + 1*x112 + 1*x210 + 1*x211 + 1*x212 = 1
1*x020 + 1*x021 + 1*x022 + 1*x120 + 1*x121 + 1*x122 + 1*x220 + 1*x221 + 1*x222 = 1
1*x000 + 1*x010 + 1*x020 + 1*x100 + 1*x110 + 1*x120 + 1*x200 + 1*x210 + 1*x220 = 1
1*x001 + 1*x011 + 1*x021 + 1*x101 + 1*x111 + 1*x121 + 1*x201 + 1*x211 + 1*x221 = 1
1*x002 + 1*x012 + 1*x022 + 1*x102 + 1*x112 + 1*x122 + 1*x202 + 1*x212 + 1*x222 = 1
801*x000 + 391*x001 + 705*x002 + 596*x010 + 472*x011 + 821*x012 + 490*x020 + 261*x021 + 710*x022 + 590*x100 + 322*x101 + 373*x102 + 496*x110 + 514*x111 + 600*x112 + 492*x120 + 405*x121 + 591*x122 + 632*x200 + 535*x201 + 533*x202 + 375*x210 + 564*x211 + 597*x212 + 288*x220 + 372*x221 + 505*x222 => max";
 
            var les = LinearEquationSystem.Parse(task);
            var solver = new SimplexSolver();
            var res = solver.Solve(les);
 
            foreach (var pair in res)
                Console.WriteLine("{0} = {1:0.00}", pair.Key, pair.Value);
 
            Console.ReadLine();
        }
    }