Re: Linear optimization solution

petter — Sat, 16 May 2015 20:44:51 GMT

This might be worth having a look at:

or even this:

This last one is Simplex implemented in JavaScript and it works perfectly well in the Module Macro in QlikView if you turn on JScript:

You can copy and paste this into the Macro Module editor:

//-------------------------------------------------------------------

// SimplexJS

// https://github.com/IainNZ/SimplexJS

// A linear and mixed-integer linear programming solver.

// By Iain Dunning (c)

// Licensed under the MIT License.

//-------------------------------------------------------------------

// Place everything under the SimplexJS namespace

var SimplexJS = {

// Solver status constants

INFEASIBLE : 0,

OPTIMAL : 1,

UNBOUNDED : 2,

// Parameters

MAXITERATIONS: Infinity,

//---------------------------------------------------------------

// ModelCopy

// 'Deep copies' a model by creating a copy of all children

// (e.g. A, b, c, ...). By default, copying an ojbect just

// creates a 'shallow' copy, where the two objects would both

// refer to the same matrices.

// Variations of this code can be found everywhere with a quick

// Google search - I'm not sure exactly where I got it.

ModelCopy : function(model) {

if (model == null || typeof(model) != 'object') return model;

var newModel = new model.constructor();

for (var key in model) newModel[key] = SimplexJS.ModelCopy(model[key]);

return newModel;

//---------------------------------------------------------------

// SolveMILP

// Uses branch-and-bound to solve a (mixed-integer) linear

// program. This implementation is about as basic as it gets!

// Assumes the root model is in the form

// min c.x

// st A.x = b

// l <= x <= u

// Some x are integer, and b >= 0

SolveMILP : function(rootModel, maxNodes, log) {

if (maxNodes === undefined) maxNodes = Infinity;

if (log === undefined) log = [];

// Track the best integer solution found

var bestFeasible = Infinity;

var bestFeasibleX = new Array(rootModel.n);

// Branch on most fractional variable

var mostFracIndex, mostFracValue, fracValue;

// Create and start branch-and-bound tree

var unsolvedLPs = new Array();

rootModel.solved = false;

unsolvedLPs.push(rootModel);

var nodeCount = 0;

// Process unsolved nodes on the b-and-b tree

while (unsolvedLPs.length >= 1) {

nodeCount += 1;

model = unsolvedLPs.shift();

// Stop if we reached max node count

if (nodeCount >= maxNodes) {

log.push("Tried to solve node #"+nodeCount.toString()+" which is >= max = "+maxNodes.toString());

unsolvedLPs = [];

rootModel.status = SimplexJS.INFEASIBLE;

return;

}

// Solve the LP at this node

log.push("Solving node #"+nodeCount.toString()+", nodes on tree: "+(unsolvedLPs.length+1).toString());

SimplexJS.PrimalSimplex(model, log);

if (model.status == SimplexJS.INFEASIBLE) {

// LP is infeasible, fathom it

log.push("Node infeasible, fathoming.");

continue;

}

log.push("LP solution at node = "+model.z.toString());

//console.log(model.x);

// Is this worse than the best integer solution?

if (model.z > bestFeasible) {

// Fathom

log.push("LP solution worse than best integer feasible, fathoming.");

continue;

}

// Check integrality of LP solution

mostFracIndex = -1;

mostFracValue = 0;

for (i = 0; i < model.n; i++) {

// Is this variable integer?

if (model.xINT) {

// Is it fractional?

if (Math.abs(Math.floor(model.x) - model.x) > 0.0001) {

// Does not satisfy integrality - will need to branch

fracValue = Math.min( Math.abs(Math.floor(model.x) - model.x),

Math.abs(Math.ceil (model.x) - model.x) );

if (fracValue > mostFracValue) {

mostFracIndex = i;

mostFracValue = fracValue;

}

// Did we find any fractional ints?

if (mostFracIndex == -1) {

// No fractional ints - update best feasible?

log.push("Node is integer feasible.");

if (model.z < bestFeasible) {

log.push("Best integer feasible was "+bestFeasible.toString()+", is now = "+model.z.toString());

bestFeasible = model.z;

for (i = 0; i < model.n; i++) bestFeasibleX = model.x;

}

} else {

// Some fractional - create two new LPs to solve

log.push("Node is fractional, branching on most fractional variable, "+mostFracIndex.toString());

downBranchModel = SimplexJS.ModelCopy(model);

downBranchModel.xUB[mostFracIndex] = Math.floor(downBranchModel.x[mostFracIndex])

downBranchModel.z = model.z;

unsolvedLPs.push(downBranchModel);

upBranchModel = SimplexJS.ModelCopy(model);

upBranchModel.xLB[mostFracIndex] = Math.ceil(upBranchModel.x[mostFracIndex])

upBranchModel.z = model.z;

unsolvedLPs.push(upBranchModel);

}

// How did it go?

rootModel.nodeCount = nodeCount;

if (bestFeasible < Infinity) {

// Done!

log.push("All nodes solved or fathomed - integer solution found");

rootModel.x = bestFeasibleX;

rootModel.z = bestFeasible;

rootModel.status = SimplexJS.OPTIMAL;

} else {

log.push("All nodes solved or fathomed - NO integer solution found");

rootModel.status = SimplexJS.INFEASIBLE;

}

//---------------------------------------------------------------

// PrimalSimplex

// Uses the revised simplex method with bounds to solve the

// primal of a linear program. Assumes the model is in the form:

// min c.x

// st A.x = b

// l <= x <= u

// Some x are integer, and b >= 0

PrimalSimplex : function(model, log) {

if (log === undefined) log = [];

// Get shorter names

A=model.A; b=model.b; c=model.c;

m=model.m; n=model.n;

xLB=model.xLB; xUB=model.xUB;

// Define some temporary variables we will need for RSM

var i, j;

var varStatus = new Array(n + m);

var basicVars = new Array(m);

var Binv = new Array(m);

for (i = 0; i < m; i++) { Binv = new Array(m); }

var cBT = new Array(m);

var pi = new Array(m);

var rc = new Array(n);

var BinvAs = new Array(m);

// Some useful constants

var BASIC = 0, NONBASIC_L = +1, NONBASIC_U = -1;

var TOL = 0.000001;

// The solution

var x = new Array(n + m), z, status;

// Create the initial solution to Phase 1

// - Real variables

for (i = 0; i < n; i++) {

var absLB = Math.abs(xLB);

var absUB = Math.abs(xUB);

x = (absLB < absUB) ? xLB : xUB ;

varStatus = (absLB < absUB) ? NONBASIC_L : NONBASIC_U;

}

// - Artificial variables

for (i = 0; i < m; i++) {

x[i+n] = b;

// Some of the real variables might be non-zero, so need

// to reduce x[artificials] accordingly

for (j = 0; j < n; j++) { x[i+n] -= A * x; }

varStatus[i+n] = BASIC;

basicVars = i+n;

}

// - Basis

for (i = 0; i < m; i++) { cBT = +1.0; }

for (i = 0; i < m; i++) {

for (j = 0; j < m; j++) {

Binv = (i == j) ? 1.0 : 0.0;

}

// Being simplex iterations

var phaseOne = true, iter = 0;

while (true) {

iter++;

if (iter >= SimplexJS.MAXITERATIONS) {

log.push("PrimalSimplex: Reached MAXITERATIONS="+iter.toString());

z = 0.0;

for (i = 0; i < n; i++) z += c * x;

model.z = z; //Infinity;

model.x = x;

break;

}

//---------------------------------------------------------------------

// Step 1. Duals and reduced Costs

//console.log(Binv)

for (i = 0; i < m; i++) {

pi = 0.0;

for (j = 0; j < m; j++) {

pi += cBT * Binv

}

//console.log(pi);

for (j = 0; j < n; j++) {

rc = phaseOne ? 0.0 : c;

for (i = 0; i < m; i++) {

rc -= pi * A;

}

//console.log(rc);

//---------------------------------------------------------------------

// Step 2. Check optimality and pick entering variable

var minRC = -TOL, s = -1;

for (i = 0; i < n; i++) {

// If NONBASIC_L (= +1), rc must be negative (< 0) -> +rc < -TOL

// If NONBASIC_U (= -1), rc must be positive (> 0) -> -rc < -TOL

// -> +rc > +TOL

// If BASIC (= 0), can't use this rc -> 0 * rc < -LPG_TOL -> alway FALSE

// Then, by setting initial value of minRC to -TOL, can collapse this

// check and the check for a better RC into 1 IF statement!

if (varStatus * rc < minRC) {

minRC = varStatus * rc;

s = i;

}

//console.log(minRC, s);

// If no entering variable

if (s == -1) {

if (phaseOne) {

//console.log("Phase one optimal")

z = 0.0;

for (i = 0; i < m; i++) z += cBT * x[basicVars];

if (z > TOL) {

//console.log("Phase 1 objective: z = ", z, " > 0 -> infeasible!");

log.push("PrimalSimplex: Phase 1 objective: z = "+z.toString()+" > 0 -> infeasible!");

model.status = SimplexJS.INFEASIBLE;

break;

} else {

//log.push("Transitioning to phase 2");

phaseOne = false;

for (i = 0; i < m; i++) {

cBT = (basicVars < n) ? (c[basicVars]) : (0.0);

}

continue;

}

} else {

model.status = SimplexJS.OPTIMAL;

z = 0.0;

for (i = 0; i < n; i++) {

z += c * x;

}

model.z = z;

model.x = x;

//console.log("Optimality in Phase 2!",z);

log.push("PrimalSimplex: Optimality in Phase 2: z = "+z.toString());

//console.log(x);

break;

}

//---------------------------------------------------------------------

// Step 3. Calculate BinvAs

for (i = 0; i < m; i++) {

BinvAs = 0.0;

for (k = 0; k < m; k++) BinvAs += Binv * A;

}

//console.log(BinvAs);

//---------------------------------------------------------------------

// Step 4. Ratio test

var minRatio = Infinity, ratio = 0.0, r = -1;

var rIsEV = false;

// If EV is...

// NBL, -> rc ~~< 0 -> want to INCREASE x~~

// NBU, -> rc ~~> 0 -> want to DECREASE x~~

// Option 1: Degenerate iteration

ratio = xUB ~~- xLB;~~

if (ratio <= minRatio) { minRatio = ratio; r = -1; rIsEV = true; }

// Option 2: Basic variables leaving basis

for (i = 0; i < m; i++) {

j = basicVars;

var jLB = (j >= n) ? 0.0 : xLB;

var jUB = (j >= n) ? Infinity : xUB;

if (-1*varStatus*BinvAs > +TOL) { // NBL: BinvAs < 0, NBU: BinvAs > 0

ratio = (x - jUB) / (varStatus*BinvAs);

if (ratio <= minRatio) { minRatio = ratio; r = i; rIsEV = false; }

}

if (+1*varStatus*BinvAs > +TOL) { // NBL: BinvAs > 0, NBU: BinvAs < 0

ratio = (x - jLB) / (varStatus*BinvAs);

if (ratio <= minRatio) { minRatio = ratio; r = i; rIsEV = false; }

}

// Check ratio

if (minRatio >= Infinity) {

if (phaseOne) {

// Not sure what this means - nothing good!

//console.log("Something bad happened");

log.push("PrimalSimplex: Something bad happened in Phase 1...");

break;

} else {

// PHASE 2: Unbounded!

model.status = SimplexJS.UNBOUNDED;

//console.log("Unbounded in Phase 2!");

log.push("PrimalSimplex: Unbounded in Phase 2!");

break;

}

//---------------------------------------------------------------------

// Step 5. Update solution and basis

x ~~+= varStatus * minRatio;~~

for (i = 0; i < m; i++) x[basicVars] -= varStatus * minRatio * BinvAs;

if (!rIsEV) {

// Basis change! Update Binv, flags

// RSM tableau: [Binv B | Binv | Binv As]

// -> GJ pivot on the BinvAs column, rth row

var erBinvAs = BinvAs;

// All non-r rows

for (i = 0; i < m; i++) {

if (i != r) {

var eiBinvAsOvererBinvAs = BinvAs / erBinvAs;

for (j = 0; j < m; j++) {

Binv -= eiBinvAsOvererBinvAs * Binv

}

// rth row

for (j = 0; j < m; j++) Binv /= erBinvAs;

// Update status flags

varStatus ~~= BASIC;~~

if (basicVars < n) {

if (Math.abs(x[basicVars] - xLB[basicVars]) < TOL) varStatus[basicVars] = NONBASIC_L;

if (Math.abs(x[basicVars] - xUB[basicVars]) < TOL) varStatus[basicVars] = NONBASIC_U;

} else {

if (Math.abs(x[basicVars] - 0.00000) < TOL) varStatus[basicVars] = NONBASIC_L;

if (Math.abs(x[basicVars] - Infinity) < TOL) varStatus[basicVars] = NONBASIC_U;

}

cBT = phaseOne ? 0.0 : c;

basicVars = s;

} else {

// Degenerate iteration

if (varStatus ~~== NONBASIC_L) { varStatus ~~= NONBASIC_U; }~~~~

else { varStatus ~~= NONBASIC_L; }~~

}

//---------------------------------------------------------------------

//console.log(x);

}

};

function TestPrimalSimplex() {

var test = new Object();

test.A = [[ 2, 1, 1, 0],

[20, 1, 0, 1]];

test.b = [40, 100];

test.c = [-10, -1, 0, 0];

test.m = 2;

test.n = 4;

test.xLB = [2, 0, 0, 0];

test.xUB = [3, Infinity, Infinity, Infinity];

SimplexJS.PrimalSimplex(test);

qvlib.MsgBox(test.x, test.z);

// Should be 3, 34, 0, 6

}

function TestBandB() {

var test = new Object();

test.A = [[ 1, 1, 0, 1, 0, 0],

[ 0, 1, 1, 0, 1, 0],

[.5,.5, 1, 1, 0, 1]];

test.b = [ 1, 1, 1];

test.c = [-1,-1,-1, 0, 0, 0];

test.m = 3;

test.n = 6;

test.xLB = [0, 0, 0, 0, 0, 0];

test.xUB = [Infinity, Infinity, Infinity, Infinity, Infinity, Infinity];

test.xINT = [true, true, true, false, false, false];

SimplexJS.SolveMILP(test);

qvlib.MsgBox(test.x, test.z);

// Should be 1, 0, 0, 0, 1, 0.5, z=-1

// or 0, 1, 0, 0, 0, 0.5, z=-1

}

function BothTests() {

qvlib.MsgBox("Testing Primal");

TestPrimalSimplex();

qvlib.MsgBox("Should be 3, 34, 0, 6");

qvlib.MsgBox("Testing Branch and Bound");

TestBandB();

qvlib.MsgBox(" Should be 1, 0, 0, 0, 1, 0.5, z=-1");

qvlib.MsgBox(" or 0, 1, 0, 0, 0, 0.5, z=-1");

}

topic Linear optimization solution in QlikView

Linear optimization solution

Re: Linear optimization solution

Re: Linear optimization solution

Re: Linear optimization solution