Trying to solve a linear program with high precision requirements and am interested in both the primal and the dual solution. Prefer solving the primal problem for convenience and solving the dual problem directly is much slower.
The problem is that beyond a low level of precision in the inputs, LinearOptimization produces a dual solution that is not feasible and not close to feasible.
The smallest version of this I can find is:
Minimize error
s.t.
w[3, 1] == 1,
w[4, 1] == 1,
error - 0.711710110526 w[3, 1] + 0.889887251073 w[4, 1] >= 0.182039226889,
error - 0.988334554323 w[3, 1] + 1.16651169487 w[4, 1] >= 0.182039226889,
error + 0.711710110526 w[3, 1] - 0.889887251073 w[4, 1] >= -0.185931432889,
error + 0.988334554323 w[3, 1] - 1.16651169487 w[4, 1] >= -0.185931432889,
w[3, 1] >= 0, w[4, 1] >= 0
Code to produce that:
(* High precision coefficients I want to work with *)
frules = {f[10,3,1]->0.711710110525816048889380805190494883846397745029174914215781022415643131249209751570622867226699`77.79196577853726,
f[11,3,1]->0.988334554323414084616848317521162833676376991758444880531217652747073397938602080305120566993438`77.9550710812633,
f[10,4,1]->0.889887251073471354449280136393810439744021905005001130029585130019076635434684301260170902240083`78.52043561129913,
f[11,4,1]->1.16651169487106939017674764872447838957400115173427109634502176035050690212407662999466860200682`78.73394921448386};
kays ={k[3,4]->0.0015630489999999999488344837317299607093445956707000732421875`80.,
k[4,3]->0.00232915700000000010227996227740732138045132160186767578125`80.};
sol = {F[3]->0.162544067765155662247056777958811870662084298098888628443782481342830102638614712905757204334513`78.62152715754068,
F[4]->-0.021058208123808663871846254592786671674207332500808453607955220577185781639628248477835432611615`78.22963096482684};
(* Make the problem as simple as possible *)
npaths[3]=1;
npaths[4]=1;
(* Trim the data to the minimum acceptable precision *)
frules = (frules/.Rule[exp_,val_]:>Rule[exp,N[val,12]]);
(* The dual variables for mapping solutions *)
dualys = Flatten[Table[With[{modei=modei,modej=modej, pti=pti},
If[modei==modej, {}, y[modei,modej,pti]]], (* k-constraints for each i, j and pt *)
{modei, 3, 4}, {modej, 3, 4}, {pti, 10, 11}]];
dualpis = Table[With[{modei=modei}, pi[modei]],{modei, 3, 4}];
(* Solve the primal problem *)
{val, {ys, pis}, ws} = LinearOptimization[error, (* Objective *)
(* Constraints *)
Join[Table[With[{modei=modei},
\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(pathi = 1\), \(npaths[modei]\)]\(w[modei, pathi]\)\)==1],(* convex combination of paths for each mode *)
{modei, 3, 4}],
Flatten[Table[With[{modei=modei,modej=modej, pti=pti},
If[modei==modej, {},
(error - \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(pathi = 1\), \(npaths[modei]\)]\(w[modei, pathi] f[pti, modei, pathi]\)\)
+ \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(pathj = 1\), \(npaths[modej]\)]\(w[modej, pathj] f[pti, modej, pathj]\)\)/.frules)
>= N[F[modei]-F[modej]-k[modei,modej]/.sol/.kays,12]]], (* k-constraints for each i, j and pt *)
{modei, 3, 4}, {modej, 3, 4}, {pti, 10, 11}]],
Flatten[Table[With[{modei=modei,pathi=pathi}, w[modei,pathi]>=0],{modei, 3, 4},{pathi, 1, npaths[modei]}]]], (* We are only working on the paths from s[i,\[Tau][i]] to \[CapitalTheta] here *)
(* Variables *)
Join[{error},Flatten[Table[With[{modei=modei,pathi=pathi}, w[modei,pathi]],{modei, 3, 4},{pathi, 1, npaths[modei]}]]],
{"PrimalMinimumValue","DualMaximizer", "PrimalMinimizerRules"}];
(* Map the dual variables to their values *)
dualsol = Flatten[Join[Thread[dualys->ys[[1;;Length[dualys]]]],
Thread[dualpis->pis]]];
(* Share the results *)
Column[Join[{Row[{Style["Primal Solution",Bold]}]},
ws,
{Row[{Style["Dual Solution",Bold]}]},
dualsol]]
Output from LinearOptimization:
Primal Solution
error -> 0.003862086,
w[3, 1] -> 1.0000000,
w[4, 1] -> 1.0000000,
Dual Solution,
y[3, 4, 10] -> 1.0000000,
y[3, 4, 11] -> 0,
y[4, 3, 10] -> 0,
y[4, 3, 11] -> 0,
pi[3] -> 0,
pi[4] -> 0.8898873
The dual problem (in Minimize and >= form):
Minimize -pi[3] - pi[4] - 0.182039226889 y[3, 4, 10] - 0.182039226889 y[3, 4, 11] + 0.185931432889 y[4, 3, 10] + 0.185931432889 y[4, 3, 11],
s.t.
y[3, 4, 10] + y[3, 4, 11] + y[4, 3, 10] + y[4, 3, 11] == 1,
-pi[3] - 0.711710110526 (-y[3, 4, 10] + y[4, 3, 10]) - 0.988334554323 (-y[3, 4, 11] + y[4, 3, 11]) >= 0,
-pi[4] - 0.889887251073 (y[3, 4, 10] - y[4, 3, 10]) - 1.16651169487 (y[3, 4, 11] - y[4, 3, 11]) >= 0,
y[3, 4, 10] >= 0, y[3, 4, 11] >= 0, y[4, 3, 10] >= 0, y[4, 3, 11] >= 0
Code to solve that:
dualout = LinearOptimization[-\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(modei = 3\), \(4\)]\(pi[modei]\)\)- \!\(
\*UnderoverscriptBox[\(\[Sum]\), \(pti = 10\), \(11\)]\(
\*UnderoverscriptBox[\(\[Sum]\), \(modei = 3\), \(4\)]\(
\*UnderoverscriptBox[\(\[Sum]\), \(modej = 3\), \(4\)]N[\(\((F[modei] - F[modej] - k[modei, modej])\) /. sol\) /. kays, 12] If[modei == modej, 0, y[modei, modej, pti]]\)\)\),
Join[
{\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(pti = 10\), \(11\)]\(
\*UnderoverscriptBox[\(\[Sum]\), \(modei = 3\), \(4\)]\(
\*UnderoverscriptBox[\(\[Sum]\), \(modej = 3\), \(4\)]If[modei == modej, \ 0, y[modei, modej, pti]]\)\)\) ==1},
Flatten[Table[With[{modei=modei},
-pi[modei]-\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(pti = 10\), \(11\)]\(
\*UnderoverscriptBox[\(\[Sum]\), \(modej = 3\), \(4\)]\((y[modej, modei, pti] - y[modei, modej, pti])\) f[pti, modei, 1]\)\)>=0/.frules],{modei,3, 4}]],
Flatten[Table[With[{modei=modei,modej=modej,pti=pti}, If[modei==modej, {}, y[modei,modej,pti]>= 0]],{modei, 3, 4}, {modej, 3, 4}, {pti, 10, 11}]]],
Join[Table[With[{modei=modei}, pi[modei]], {modei, 3, 4}], Flatten[Table[With[{modei=modei,modej=modej, pti=pti}, If[modei==modej,{},y[modei,modej,pti]]],{modei, 3, 4}, {modej, 3, 4}, {pti, 10, 11}]]],
{"PrimalMinimumValue", "PrimalMinimizerRules","DualMaximizer"},
Method -> "RevisedSimplex"
];
Column[Join[{Row[{Style["Dual Solution",Bold]}]},
Flatten[dualout[[1;;2]]]]]
The results from solving the dual directly:
Dual Solution
objective value: -0.00386209
pi[3] -> 0.71171011
pi[4] -> -0.88988725
y[3, 4, 10] -> 1.00000000
y[3, 4, 11] -> 0
y[4, 3, 10] -> 0
y[4, 3, 11] -> 0
Is the dual solution produced when solving the primal problem feasible? No! In particular, the third constraint:
-pi[4]-0.889887251073 (y[3,4,10]-y[4,3,10])-1.16651169487 (y[3,4,11]-y[4,3,11])>=0
is not close to being satisfied:
-pi[4]-f[10,4,1] (y[3,4,10]-y[4,3,10])-f[11,4,1] (y[3,4,11]-y[4,3,11])/.frules/.dualsol
yields -1.779775.
Including WorkingPrecision->12 as an option in solving the primal problem does not change this. Reducing the precision of the formulation to, for example, 6 digits resolves the issue, but is too crude for the purposes at hand. Am I mistaken in believing that LinearOptimization can support arbitrary precision as suggested by the WorkingPrecision option? Indeed, converting parameters to exact numbers and increasing WorkingPrecision does not change this. Suggestions?
