algorithm - 手动编写线性编程练习

Question

我一直在课堂上通过绘制线性规划问题来解决它们，但我想知道如何为特定问题编写程序来为我解决它。如果变量或约束太多，我永远无法通过绘图来做到这一点。

示例问题，在约束条件下最大化 5x + 3y：

• 5x - 2y >= 0
• x + y <= 7
• x <= 5
• x >= 0
• y >= 0

我绘制了这个图，得到了一个有 3 个角的可见区域。x=5 y=2 是最佳点。

我如何把它变成代码？我知道单纯形法。并且非常重要的是，所有 LP 问题都会被编码在同一个结构中吗？蛮力会起作用吗？

Answer 1

如果您搜索，您会发现相当多的 Simplex 实现。

除了评论中提到的（C 中的数字食谱），您还可以找到：

1. Google 自己的 Simplex-Solver
2. 然后是硬币或
3. GNU 有自己的GLPK
4. 如果你想要一个 C++ 实现，Google Code中的这个实际上是可以访问的。
5. R 中有许多实现，包括引导包。（在 R 中，您可以通过键入不带括号的函数来查看函数的实现。）

要解决您的其他两个问题：

1. 所有 LP 的编码方式都相同吗？是的，可以编写一个通用的 LP 求解器来加载和求解任何 LP。（有读取 LP 之类的行业标准格式mps和.lp

2. 蛮力会起作用吗？请记住，许多公司和大型组织花费很长时间来微调求解器。有些 LP 具有许多求解器会尝试利用的有趣特性。此外，某些计算可以并行解决。该算法是指数的，因此在大量变量/约束下，蛮力将不起作用。

希望有帮助。

Answer 2

我昨天写的是 matlab，如果你使用 Eigen 库或使用 std::vector 的 std::vector 编写自己的矩阵类，它可以很容易地转录为 C++

function [x, fval] = mySimplex(fun, A, B, lb, up)

%Examples paramters to show that the function actually works 

% sample set 1 (works for this data set)

% fun = [8 10 7];
% A = [1 3 2; 1 5 1];
% B = [10; 8];
% lb = [0; 0; 0];
% ub = [inf; inf; inf];

% sample set 2 (works for this data set)

fun = [7 8 10];
A = [2 3 2; 1 1 2];
B = [1000; 800];
lb = [0; 0; 0];
ub = [inf; inf; inf];


% generate a new slack variable for every row of A 

numSlackVars = size(A,1); % need a new slack variables for every row of A 

% Set up tableau to store algorithm data 
tableau = [A; -fun];

tableau = [tableau, eye(numSlackVars + 1)];

lastCol = [B;0];

tableau = [tableau, lastCol];

% for convienience sake, assign the following: 

numRows = size(tableau,1);
numCols = size(tableau,2);

% do simplex algorithm 

% step 0: find num of negative entries in bottom row of tableau 

numNeg = 0; % the number of negative entries in bottom row

for i=1:numCols 
    if(tableau(numRows,i) < 0)
        numNeg = numNeg + 1;
    end
end

% Remark: the number of negatives is exactly the number of iterations needed in the
% simplex algorithm 

for iterations = 1:numNeg 
    % step 1: find minimum value in last row 
    minVal = 10000; % some big number 
    minCol = 1; % start by assuming min value is the first element 
    for i=1:numCols
        if(tableau(numRows, i) < minVal)
            minVal = tableau(size(tableau,1), i);
            minCol = i; % update the index corresponding to the min element 
        end
    end 

    % step 2: Find corresponding ratio vector in pivot column 
    vectorRatio = zeros(numRows -1, 1);
    for i=1:(numRows-1) % the size of ratio vector is numCols - 1
        vectorRatio(i, 1) = tableau(i, numCols) ./ tableau(i, minCol);
    end 

    % step 3: Determine pivot element by finding minimum element in vector
    % ratio

    minVal = 10000; % some big number 
    minRatio = 1; % holds the element with the minimum ratio 

    for i=1:numRows-1
        if(vectorRatio(i,1) < minVal)
            minVal = vectorRatio(i,1);
            minRatio = i;
        end 
    end 

    % step 4: assign pivot element 

    pivotElement = tableau(minRatio, minCol);

    % step 5: perform pivot operation on tableau around the pivot element 

    tableau(minRatio, :) = tableau(minRatio, :) * (1/pivotElement);

    % step 6: perform pivot operation on rows (not including last row)

    for i=1:size(vectorRatio,1)+1 % do last row last 
        if(i ~= minRatio) % we skip over the minRatio'th element of the tableau here 
            tableau(i, :) = -tableau(i,minCol)*tableau(minRatio, :) +  tableau(i,:);
        end
    end
end 

% Now we can interpret the algo tableau 

numVars = size(A,2); % the number of cols of A is the number of variables 

x = zeros(size(size(tableau,1), 1)); % for efficiency 

% Check for basicity 
for col=1:numVars
    count_zero = 0;
    count_one = 0;
    for row = 1:size(tableau,1)
        if(tableau(row,col) < 1e-2)
            count_zero = count_zero + 1;
        elseif(tableau(row,col) - 1 < 1e-2)
            count_one = count_one + 1;
            stored_row = row; % we store this (like in memory) column for later use 
        end
    end
    if(count_zero == (size(tableau,1) -1) && count_one == 1) % this is the case where it is basic 
        x(col,1) = tableau(stored_row, numCols);
    else 
        x(col,1) = 0; % this is the base where it is not basic 
    end
end

% find function optimal value at optimal solution 
fval = x(1,1) * fun(1,1); % just needed for logic to work here 
for i=2:numVars 
    fval = fval + x(i,1) * fun(1,i);
end


end