MDP solution via linear programming for explicit engine.

prism-tests/functionality/verify/mdps/reach/mdp_simple.nm.props.args

@@ -5,6 +5,7 @@
 -ex -gs
 -ex -politer
 -ex -modpoliter
+-ex -lp
 -s -ii
 -m -ii
 -h -ii

prism-tests/functionality/verify/mdps/rewards/rewpoliter.nm.props.args

@@ -4,3 +4,4 @@
 -valiter -h
 -valiter -ex
 -gs -ex
+-gs -lp

prism-tests/pmc/lec13and14mdp.nm

@@ -16,3 +16,8 @@ s:[0..3];
 endmodule
 
 label "a" = s=2;
+
+rewards
+	[] true : 0.6;
+	true : 0.4;
+endrewards

prism-tests/pmc/lec13and14mdp.nm.props

@@ -12,3 +12,8 @@
 	filter(forall, P<0.1 [ F "a" ] <=> false);
 	// RESULT: true
 	filter(forall, P>0 [ F "a" ] <=> (s=0|s=1|s=2));
+
+	// RESULT: 5/3
+	Rmin=? [ F "a" ];
+	// RESULT: Infinity
+	Rmax=? [ F "a" ];

prism/src/explicit/MDPModelChecker.java

@@ -49,8 +49,12 @@
 import explicit.rewards.MDPRewards;
 import explicit.rewards.Rewards;
 import io.ModelExportFormat;
+import lpsolve.LpSolve;
+import lpsolve.LpSolveException;
 import parser.ast.Expression;
 import parser.type.TypeDouble;
+import prism.Accuracy;
+import prism.Accuracy.AccuracyLevel;
 import prism.AccuracyFactory;
 import prism.OptionsIntervalIteration;
 import prism.PrismComponent;
@@ -318,12 +322,6 @@ public ModelCheckerResult computeReachProbs(MDP<Double> mdp, BitSet remain, BitS
 
 		boolean doPmaxQuotient = this.doPmaxQuotient;
 
-		// Switch to a supported method, if necessary
-		if (mdpSolnMethod == MDPSolnMethod.LINEAR_PROGRAMMING) {
-			mdpSolnMethod = MDPSolnMethod.GAUSS_SEIDEL;
-			mainLog.printWarning("Switching to MDP solution method \"" + mdpSolnMethod.fullName() + "\"");
-		}
-
 		// Check for some unsupported combinations
 		if (mdpSolnMethod == MDPSolnMethod.VALUE_ITERATION && valIterDir == ValIterDir.ABOVE) {
 			if (!(precomp && prob0))
@@ -353,6 +351,14 @@ public ModelCheckerResult computeReachProbs(MDP<Double> mdp, BitSet remain, BitS
 				throw new PrismException("Policy iteration methods cannot be passed 'known' values for some states");
 			}
 		}
+		if (mdpSolnMethod == MDPSolnMethod.LINEAR_PROGRAMMING) {
+			if (!(precomp && prob0)) {
+				throw new PrismNotSupportedException("Prob0 precomputation must be enabled for linear programming");
+			}
+			if (!min && genStrat) {
+				throw new PrismNotSupportedException("Currently, explicit engine does not support strategy generation for linear programming and Pmax");
+			}
+		}
 
 		if (doPmaxQuotient && min) {
 			// for Pmin, don't do quotient
@@ -535,6 +541,9 @@ protected ModelCheckerResult computeReachProbsNumeric(MDP<Double> mdp, MDPSolnMe
 			}
 			res = computeReachProbsModPolIter(mdp, no, yes, min, strat);
 			break;
+		case LINEAR_PROGRAMMING:
+			res = computeReachProbsLP(mdp, no, yes, min, strat);
+			break;
 		default:
 			throw new PrismException("Unknown MDP solution method " + mdpSolnMethod.fullName());
 		}
@@ -1189,6 +1198,124 @@ protected ModelCheckerResult computeReachProbsModPolIter(MDP<Double> mdp, BitSet
 		return res;
 	}
 
+	/**
+	 * Compute reachability probabilities using linear programming.
+	 * @param mdp: The MDP
+	 * @param no: Probability 0 states
+	 * @param yes: Probability 1 states
+	 * @param min: Min or max probabilities (true=min, false=max)
+	 * @param strat Storage for (memoryless) strategy choice indices (ignored if null)
+	 */
+	protected ModelCheckerResult computeReachProbsLP(MDP<Double> mdp, BitSet no, BitSet yes, boolean min, int[] strat) throws PrismException
+	{
+		double[] soln = null;
+		long timer;
+
+		// Start solution
+		timer = System.currentTimeMillis();
+		mainLog.println("Starting linear programming (" + (min ? "min" : "max") + ")...");
+
+		// Store num states
+		int n = mdp.getNumStates();
+
+		// Determine set of states actually need to perform computation for
+		BitSet unknown = new BitSet();
+		unknown.set(0, n);
+		unknown.andNot(yes);
+		unknown.andNot(no);
+
+		try {
+			// Initialise LP solver
+			LpSolve solver = LpSolve.makeLp(0, n);
+			try {
+				solver.setVerbose(lpsolve.LpSolve.CRITICAL);
+				solver.setAddRowmode(true);
+				// Set up arrays for passing LP to solver
+				double[] row = new double[n + 1];
+				int[] colno = new int[n + 1];
+				// Set objective function
+				if (min) {
+					solver.setMaxim();
+				} else {
+					solver.setMinim();
+				}
+				for (int s : new IterableBitSet(unknown)) {
+					row[s + 1] = 1.0;
+				}
+				solver.setObjFn(row);
+				// Add constraints
+				for (int s = 0; s < n; s++) {
+					solver.setBounds(s + 1, 0, 1);
+					if (yes.get(s)) {
+						row[0] = 1.0;
+						colno[0] = s + 1;
+						solver.addConstraintex(1, row, colno, LpSolve.EQ, 1.0);
+					} else if (no.get(s)) {
+						row[0] = 1.0;
+						colno[0] = s + 1;
+						solver.addConstraintex(1, row, colno, LpSolve.EQ, 0.0);
+					} else {
+						int numChoices = mdp.getNumChoices(s);
+						for (int i = 0; i < numChoices; i++) {
+							int count = 0;
+							row[count] = 1.0;
+							colno[count] = s + 1;
+							count++;
+							Iterator<Map.Entry<Integer, Double>> iter = mdp.getTransitionsIterator(s, i);
+							while (iter.hasNext()) {
+								Map.Entry<Integer, Double> e = iter.next();
+								int t = e.getKey();
+								double p = e.getValue();
+								if (t == s) {
+									row[0] -= p;
+								} else {
+									row[count] = -p;
+									colno[count] = t + 1;
+									count++;
+								}
+							}
+							solver.addConstraintex(count, row, colno, min ? LpSolve.LE : LpSolve.GE, 0.0);
+						}
+					}
+				}
+				// Solve LP
+				solver.setAddRowmode(false);
+				//solver.printLp();
+				int lpRes = solver.solve();
+				if (lpRes == lpsolve.LpSolve.OPTIMAL) {
+					soln = solver.getPtrVariables();
+				} else {
+					String detail = lpRes == LpSolve.INFEASIBLE ? " (infeasible)" : lpRes == LpSolve.UNBOUNDED ? " (unbounded)" : "";
+					throw new PrismException("Error solving LP" + detail);
+				}
+			} finally {
+				solver.deleteLp();
+			}
+		} catch (LpSolveException e) {
+			throw new PrismException("Error solving LP: " + e.getMessage());
+		}
+
+		// We do an extra iteration of VI (using GS since we have just one solution vector)
+		// If strategy generation was requested (strat != null), this takes care of it
+		// (NB: that only works reliably for min - max needs more work - but this is disallowed earlier)
+		// It also has the side-effect of sometimes fixing small round-off errors from LP
+		mdp.mvMultGSMinMax(soln, min, unknown, false, true, strat);
+
+		// Finished solution
+		timer = System.currentTimeMillis() - timer;
+		mainLog.print("Linear programming");
+		mainLog.println(" took " + timer / 1000.0 + " seconds.");
+
+		// Return results
+		// (Note we don't add the strategy - the one passed in is already there
+		// and might have some existing choices stored for other states).
+		ModelCheckerResult res = new ModelCheckerResult();
+		res.soln = soln;
+		res.accuracy = new Accuracy(AccuracyLevel.EXACT_FLOATING_POINT);
+		res.timeTaken = timer / 1000.0;
+		return res;
+	}
+
 	/**
 	 * Construct strategy information for min/max reachability probabilities.
 	 * (More precisely, list of indices of choices resulting in min/max.)
@@ -2078,7 +2205,7 @@ public ModelCheckerResult computeReachRewards(MDP<Double> mdp, MDPRewards<Double
 		MDPSolnMethod mdpSolnMethod = this.mdpSolnMethod;
 
 		// Switch to a supported method, if necessary
-		if (!(mdpSolnMethod == MDPSolnMethod.VALUE_ITERATION || mdpSolnMethod == MDPSolnMethod.GAUSS_SEIDEL || mdpSolnMethod == MDPSolnMethod.POLICY_ITERATION)) {
+		if (!(mdpSolnMethod == MDPSolnMethod.VALUE_ITERATION || mdpSolnMethod == MDPSolnMethod.GAUSS_SEIDEL || mdpSolnMethod == MDPSolnMethod.POLICY_ITERATION || mdpSolnMethod == MDPSolnMethod.LINEAR_PROGRAMMING)) {
 			mdpSolnMethod = MDPSolnMethod.GAUSS_SEIDEL;
 			mainLog.printWarning("Switching to MDP solution method \"" + mdpSolnMethod.fullName() + "\"");
 		}
@@ -2094,6 +2221,14 @@ public ModelCheckerResult computeReachRewards(MDP<Double> mdp, MDPRewards<Double
 				throw new PrismNotSupportedException("Currently, explicit engine only supports interval iteration with value iteration or Gauss-Seidel for MDPs");
 			}
 		}
+		if (mdpSolnMethod == MDPSolnMethod.LINEAR_PROGRAMMING) {
+			if (!(precomp && prob1)) {
+				throw new PrismNotSupportedException("Prob1 precomputation must be enabled for linear programming");
+			}
+			if (!min && genStrat) {
+				throw new PrismNotSupportedException("Currently, explicit engine does not support strategy generation for linear programming and Rmax");
+			}
+		}
 
 		// Start expected reachability
 		timer = System.currentTimeMillis();
@@ -2245,6 +2380,9 @@ protected ModelCheckerResult computeReachRewardsNumeric(MDP<Double> mdp, MDPRewa
 			}
 			res = computeReachRewardsPolIter(mdp, mdpRewards, target, inf, min, strat);
 			break;
+		case LINEAR_PROGRAMMING:
+			res = computeReachRewardsLP(mdp, mdpRewards, target, inf, min, strat);
+			break;
 		default:
 			throw new PrismException("Unknown MDP solution method " + method.fullName());
 		}
@@ -2628,6 +2766,133 @@ protected ModelCheckerResult computeReachRewardsPolIter(MDP<Double> mdp, MDPRewa
 		return res;
 	}
 
+	/**
+	 * Compute expected reachability rewards using linear programming.
+	 * @param mdp: The MDP
+	 * @param mdpRewards The rewards
+	 * @param target Target states
+	 * @param inf States for which reward is infinite
+	 * @param min: Min or max rewards (true=min, false=max)
+	 * @param strat Storage for (memoryless) strategy choice indices (ignored if null)
+	 */
+	protected ModelCheckerResult computeReachRewardsLP(MDP<Double> mdp, MDPRewards<Double> mdpRewards, BitSet target, BitSet inf, boolean min, int[] strat) throws PrismException
+	{
+		double[] soln = null;
+		long timer;
+
+		// Start solution
+		timer = System.currentTimeMillis();
+		mainLog.println("Starting linear programming (" + (min ? "min" : "max") + ")...");
+
+		// Store num states
+		int n = mdp.getNumStates();
+
+		// Determine set of states actually need to perform computation for
+		BitSet unknown = new BitSet();
+		unknown.set(0, n);
+		unknown.andNot(target);
+		unknown.andNot(inf);
+
+		try {
+			// Initialise LP solver
+			LpSolve solver = LpSolve.makeLp(0, n);
+			try {
+				solver.setVerbose(lpsolve.LpSolve.CRITICAL);
+				solver.setAddRowmode(true);
+				// Set up arrays for passing LP to solver
+				double[] row = new double[n + 1];
+				int[] colno = new int[n + 1];
+				// Set objective function
+				if (min) {
+					solver.setMaxim();
+				} else {
+					solver.setMinim();
+				}
+				for (int s : new IterableBitSet(unknown)) {
+					row[s + 1] = 1.0;
+				}
+				solver.setObjFn(row);
+				// Add constraints
+				for (int s = 0; s < n; s++) {
+					solver.setBounds(s + 1, 0, Double.POSITIVE_INFINITY);
+					if (!unknown.get(s)) {
+						row[0] = 1.0;
+						colno[0] = s + 1;
+						solver.addConstraintex(1, row, colno, LpSolve.EQ, 0.0);
+					} else {
+						int numChoices = mdp.getNumChoices(s);
+						for (int i = 0; i < numChoices; i++) {
+							int count = 0;
+							row[count] = 1.0;
+							colno[count] = s + 1;
+							count++;
+							boolean toInf = false;
+							Iterator<Map.Entry<Integer, Double>> iter = mdp.getTransitionsIterator(s, i);
+							while (iter.hasNext()) {
+								Map.Entry<Integer, Double> e = iter.next();
+								int t = e.getKey();
+								double p = e.getValue();
+								if (t == s) {
+									row[0] -= p;
+								} else {
+									row[count] = -p;
+									colno[count] = t + 1;
+									count++;
+								}
+								toInf = toInf || inf.get(t);
+							}
+							// We ignore choices that may lead to an inf states
+							// (inf states are not handled by the LP so given value 0)
+							if (!toInf) {
+								double rew = mdpRewards.getStateReward(s) + mdpRewards.getTransitionReward(s, i);
+								solver.addConstraintex(count, row, colno, min ? LpSolve.LE : LpSolve.GE, rew);
+							}
+						}
+					}
+				}
+				// Solve LP
+				solver.setAddRowmode(false);
+				//solver.printLp();
+				int lpRes = solver.solve();
+				if (lpRes == lpsolve.LpSolve.OPTIMAL) {
+					soln = solver.getPtrVariables();
+				} else {
+					String detail = lpRes == LpSolve.INFEASIBLE ? " (infeasible)" : lpRes == LpSolve.UNBOUNDED ? " (unbounded)" : "";
+					throw new PrismException("Error solving LP" + detail);
+				}
+			} finally {
+				solver.deleteLp();
+			}
+		} catch (LpSolveException e) {
+			throw new PrismException("Error solving LP: " + e.getMessage());
+		}
+
+		// Set value for states with infinite rewards
+		for (int s : new IterableBitSet(inf)) {
+			soln[s] = Double.POSITIVE_INFINITY;
+		}
+
+		// We do an extra iteration of VI (using GS since we have just one solution vector)
+		// If strategy generation was requested (strat != null), this takes care of it
+		// (NB: that only works reliably for min - max needs more work - but this is disallowed earlier)
+		// It also has the side-effect of sometimes fixing small round-off errors from LP
+		mdp.mvMultRewGSMinMax(soln, mdpRewards, min, unknown, false, true, strat);
+
+		// Finished solution
+		timer = System.currentTimeMillis() - timer;
+		mainLog.print("Linear programming");
+		mainLog.println(" took " + timer / 1000.0 + " seconds.");
+
+		// Return results
+		// (Note we don't add the strategy - the one passed in is already there
+		// and might have some existing choices stored for other states).
+		ModelCheckerResult res = new ModelCheckerResult();
+		res.soln = soln;
+		res.accuracy = new Accuracy(AccuracyLevel.EXACT_FLOATING_POINT);
+		res.timeTaken = timer / 1000.0;
+		return res;
+	}
+
 	/**
 	 * Construct strategy information for min/max expected reachability.
 	 * (More precisely, list of indices of choices resulting in min/max.)