diff --git a/prism/include/DoubleVector.h b/prism/include/DoubleVector.h
index 3d7ea5a1..7c62f701 100644
--- a/prism/include/DoubleVector.h
+++ b/prism/include/DoubleVector.h
@@ -151,6 +151,14 @@ JNIEXPORT jdouble JNICALL Java_dv_DoubleVector_DV_1MinOverBDD
 JNIEXPORT jdouble JNICALL Java_dv_DoubleVector_DV_1MaxOverBDD
   (JNIEnv *, jobject, jlong, jlong, jlong, jint, jlong);
 
+/*
+ * Class:     dv_DoubleVector
+ * Method:    DV_MaxFiniteOverBDD
+ * Signature: (JJJIJ)D
+ */
+JNIEXPORT jdouble JNICALL Java_dv_DoubleVector_DV_1MaxFiniteOverBDD
+  (JNIEnv *, jobject, jlong, jlong, jlong, jint, jlong);
+
 /*
  * Class:     dv_DoubleVector
  * Method:    DV_SumOverBDD
diff --git a/prism/include/JDD.h b/prism/include/JDD.h
index 064931f7..56ba0654 100644
--- a/prism/include/JDD.h
+++ b/prism/include/JDD.h
@@ -455,6 +455,14 @@ JNIEXPORT jdouble JNICALL Java_jdd_JDD_DD_1FindMinPositive
 JNIEXPORT jdouble JNICALL Java_jdd_JDD_DD_1FindMax
   (JNIEnv *, jclass, jlong);
 
+/*
+ * Class:     jdd_JDD
+ * Method:    DD_FindMaxFinite
+ * Signature: (J)D
+ */
+JNIEXPORT jdouble JNICALL Java_jdd_JDD_DD_1FindMaxFinite
+  (JNIEnv *, jclass, jlong);
+
 /*
  * Class:     jdd_JDD
  * Method:    DD_RestrictToFirst
diff --git a/prism/include/dd_term.h b/prism/include/dd_term.h
index c627ec3a..44249c99 100644
--- a/prism/include/dd_term.h
+++ b/prism/include/dd_term.h
@@ -43,6 +43,7 @@ bool DD_EqualSupNormRel(DdManager *ddman, DdNode *dd1, DdNode *dd2, double epsil
 double DD_FindMin(DdManager *ddman, DdNode *dd);
 double DD_FindMinPositive(DdManager *ddman, DdNode *dd);
 double DD_FindMax(DdManager *ddman, DdNode *dd);
+double DD_FindMaxFinite(DdManager *ddman, DdNode *dd);
 DdNode *DD_RestrictToFirst(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars);
 bool DD_IsZeroOneMTBDD(DdManager *ddman, DdNode *dd);
 
diff --git a/prism/include/dv.h b/prism/include/dv.h
index 2b39ea1a..82584217 100644
--- a/prism/include/dv.h
+++ b/prism/include/dv.h
@@ -77,6 +77,7 @@ EXPORT void max_double_vector_mtbdd(DdManager *ddman, double *vec, DdNode *vec2,
 EXPORT double get_first_from_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd);
 EXPORT double min_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd);
 EXPORT double max_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd);
+EXPORT double max_finite_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd);
 EXPORT double sum_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd);
 EXPORT double sum_double_vector_over_mtbdd(DdManager *ddman, double *vec, DdNode *mult, DdNode **vars, int num_vars, ODDNode *odd);
 EXPORT void sum_double_vector_over_dd_vars(DdManager *ddman, double *vec, double *vec2, DdNode **vars, int num_vars, int first_var, int last_var, ODDNode *odd, ODDNode *odd2);
diff --git a/prism/src/dd/dd_term.cc b/prism/src/dd/dd_term.cc
index 24d3ec44..0ee4e38e 100644
--- a/prism/src/dd/dd_term.cc
+++ b/prism/src/dd/dd_term.cc
@@ -277,6 +277,33 @@ DdNode *dd
 
 //------------------------------------------------------------------------------
 
+// Find maximal finite terminal node (constant). Returns -infinity if there is none
+double DD_FindMaxFinite
+(
+DdManager *ddman,
+DdNode *dd
+)
+{
+	DdGen *gen;
+	DdNode *node;
+	bool rv = true;
+
+	double max_v = -std::numeric_limits<double>::infinity();
+
+	Cudd_ForeachNode(ddman, dd, gen, node) {
+		if (Cudd_IsConstant(node)) {
+			double v = Cudd_V(node);
+			if (v < std::numeric_limits<double>::infinity() && v > max_v) {
+				max_v = v;
+			}
+		}
+	}
+
+	return max_v;
+}
+
+//------------------------------------------------------------------------------
+
 double DD_FindMax
 (
 DdManager *ddman,
diff --git a/prism/src/dv/DoubleVector.cc b/prism/src/dv/DoubleVector.cc
index a861b337..5be7c6e8 100644
--- a/prism/src/dv/DoubleVector.cc
+++ b/prism/src/dv/DoubleVector.cc
@@ -398,6 +398,28 @@ jlong __jlongpointer odd
 
 //------------------------------------------------------------------------------
 
+JNIEXPORT jdouble JNICALL Java_dv_DoubleVector_DV_1MaxFiniteOverBDD
+(
+JNIEnv *env,
+jobject obj,
+jlong __jlongpointer vector,
+jlong __jlongpointer filter,
+jlong __jlongpointer vars,
+jint num_vars,
+jlong __jlongpointer odd
+)
+{
+	return (jdouble)max_finite_double_vector_over_bdd(
+		ddman,
+		jlong_to_double(vector),
+		jlong_to_DdNode(filter),
+		jlong_to_DdNode_array(vars), num_vars,
+		jlong_to_ODDNode(odd)
+	);
+}
+
+//------------------------------------------------------------------------------
+
 JNIEXPORT jdouble JNICALL Java_dv_DoubleVector_DV_1SumOverBDD
 (
 JNIEnv *env,
diff --git a/prism/src/dv/DoubleVector.java b/prism/src/dv/DoubleVector.java
index 1b66e4a2..40aa0ba3 100644
--- a/prism/src/dv/DoubleVector.java
+++ b/prism/src/dv/DoubleVector.java
@@ -265,7 +265,13 @@ public class DoubleVector
 	{
 		return DV_MaxOverBDD(v, filter.ptr(), vars.array(), vars.n(), odd.ptr());
 	}
-	
+
+	private native double DV_MaxFiniteOverBDD(long v, long filter, long vars, int num_vars, long odd);
+	public double maxFiniteOverBDD(JDDNode filter, JDDVars vars, ODDNode odd)
+	{
+		return DV_MaxOverBDD(v, filter.ptr(), vars.array(), vars.n(), odd.ptr());
+	}
+
 	// sum elements of vector according to a bdd (used for csl steady state operator)
 	private native double DV_SumOverBDD(long v, long filter, long vars, int num_vars, long odd);
 	public double sumOverBDD(JDDNode filter, JDDVars vars, ODDNode odd)
diff --git a/prism/src/dv/dv.cc b/prism/src/dv/dv.cc
index 53f513bd..315dc4e3 100644
--- a/prism/src/dv/dv.cc
+++ b/prism/src/dv/dv.cc
@@ -29,6 +29,7 @@
 #include <math.h>
 #include <new>
 #include <cstdint>
+#include <limits>
 
 // local function prototypes
 
@@ -40,6 +41,7 @@ static void max_double_vector_mtbdd_rec(DdManager *ddman, double *vec, DdNode *v
 static double get_first_from_bdd_rec(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, int level, ODDNode *odd, int64_t o);
 static double min_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, int level, ODDNode *odd, int64_t o);
 static double max_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, int level, ODDNode *odd, int64_t o);
+static double max_finite_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, int level, ODDNode *odd, int64_t o);
 static double sum_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, int level, ODDNode *odd, int64_t o);
 static double sum_double_vector_over_mtbdd_rec(DdManager *ddman, double *vec, DdNode *mult, DdNode **vars, int num_vars, int level, ODDNode *odd, int64_t o);
 static void sum_double_vector_over_dd_vars_rec(DdManager *ddman, double *vec, double *vec2, DdNode **vars, int num_vars, int level, int first_var, int last_var, ODDNode *odd, ODDNode *odd2, int64_t o, int64_t o2);
@@ -381,6 +383,44 @@ double max_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *fil
 	}
 }
 
+
+EXPORT double max_finite_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd)
+{
+	return max_finite_double_vector_over_bdd_rec(ddman, vec, filter, vars, num_vars, 0, odd, 0);
+}
+
+double max_finite_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, int level, ODDNode *odd, int64_t o)
+{
+	double d1, d2;
+	DdNode *dd;
+
+	if (level == num_vars) {
+		if (Cudd_V(filter) > 0) {
+			double v = vec[o];
+			if (v < std::numeric_limits<double>::infinity()) {
+				return -std::numeric_limits<double>::infinity();
+			} else {
+				return v;
+			}
+		}
+		else {
+			return -std::numeric_limits<double>::infinity();
+		}
+	}
+	else {
+		d1 = d2 = -HUGE_VAL;
+		if (odd->eoff > 0) {
+			dd = (filter->index > vars[level]->index) ? filter : Cudd_E(filter);
+			d1 = max_double_vector_over_bdd_rec(ddman, vec, dd, vars, num_vars, level+1, odd->e, o);
+		}
+		if (odd->toff > 0) {
+			dd = (filter->index > vars[level]->index) ? filter : Cudd_T(filter);
+			d2 = max_double_vector_over_bdd_rec(ddman, vec, dd, vars, num_vars, level+1, odd->t, o+odd->eoff);
+		}
+		return (d1 > d2) ? d1 : d2;
+	}
+}
+
 //------------------------------------------------------------------------------
 
 // sums up the elements of a double array - but only those in the bdd passed in
diff --git a/prism/src/jdd/JDD.cc b/prism/src/jdd/JDD.cc
index 53b9122c..fc13136d 100644
--- a/prism/src/jdd/JDD.cc
+++ b/prism/src/jdd/JDD.cc
@@ -441,6 +441,13 @@ JNIEXPORT jdouble JNICALL Java_jdd_JDD_DD_1FindMax(JNIEnv *env, jclass cls, jlon
 
 //------------------------------------------------------------------------------
 
+JNIEXPORT jdouble JNICALL Java_jdd_JDD_DD_1FindMaxFinite(JNIEnv *env, jclass cls, jlong __jlongpointer dd)
+{
+	return DD_FindMaxFinite(ddman, jlong_to_DdNode(dd));
+}
+
+//------------------------------------------------------------------------------
+
 
 JNIEXPORT jlong __jlongpointer JNICALL Java_jdd_JDD_DD_1RestrictToFirst(JNIEnv *env, jclass cls, jlong __jlongpointer dd, jlong __jlongpointer vars, jint num_vars)
 {
diff --git a/prism/src/jdd/JDD.java b/prism/src/jdd/JDD.java
index 012d9f05..08a32d3e 100644
--- a/prism/src/jdd/JDD.java
+++ b/prism/src/jdd/JDD.java
@@ -91,6 +91,7 @@ public class JDD
 	private static native double DD_FindMin(long dd);
 	private static native double DD_FindMinPositive(long dd);
 	private static native double DD_FindMax(long dd);
+	private static native double DD_FindMaxFinite(long dd);
 	private static native long DD_RestrictToFirst(long dd, long vars, int num_vars);
 	private static native boolean DD_IsZeroOneMTBDD(long dd);
 	// dd_info
@@ -963,7 +964,19 @@ public class JDD
 		checkForCuddError();
 		return rv;
 	}
-	
+
+	/**
+	 * Returns maximal finite positive terminal in dd, i.e.,
+	 * If there is none, returns -infinity.
+	 * <br>[ REFS: <i>none</i>, DEREFS: <i>none</i> ]
+	 */
+	public static double FindMaxFinite(JDDNode dd)
+	{
+		double rv = DD_FindMaxFinite(dd.ptr());
+		checkForCuddError();
+		return rv;
+	}
+
 	/**
 	 * returns dd restricted to first non-zero path (cube)
 	 * <br>[ REFS: <i>result</i>, DEREFS: dd ]
diff --git a/prism/src/prism/StateValues.java b/prism/src/prism/StateValues.java
index 38e31448..d334e4b2 100644
--- a/prism/src/prism/StateValues.java
+++ b/prism/src/prism/StateValues.java
@@ -137,6 +137,15 @@ public interface StateValues extends StateVector
 	 */
 	double maxOverBDD(JDDNode filter);
 
+	/**
+	 * Get the maximum finite value of those that are in the (BDD) filter.
+	 * <br>
+	 * If the filter is empty for this vector or all values for the filter are non-finite,
+	 * returns negative infinity.
+	 * <br>[ DEREFS: <i>none</i> ]
+	 */
+	public double maxFiniteOverBDD(JDDNode filter);
+
 	/**
 	 * Get the sum of those elements that are in the (BDD) filter.
 	 * If the filter is empty for this vector, returns 0.
diff --git a/prism/src/prism/StateValuesDV.java b/prism/src/prism/StateValuesDV.java
index bcb3c5e1..7fa58f40 100644
--- a/prism/src/prism/StateValuesDV.java
+++ b/prism/src/prism/StateValuesDV.java
@@ -314,6 +314,12 @@ public class StateValuesDV implements StateValues
 		return values.maxOverBDD(filter, vars, odd);
 	}
 
+	@Override
+	public double maxFiniteOverBDD(JDDNode filter)
+	{
+		return values.maxFiniteOverBDD(filter, vars, odd);
+	}
+
 	@Override
 	public double sumOverBDD(JDDNode filter)
 	{
diff --git a/prism/src/prism/StateValuesMTBDD.java b/prism/src/prism/StateValuesMTBDD.java
index 100cdd40..4c99e39d 100644
--- a/prism/src/prism/StateValuesMTBDD.java
+++ b/prism/src/prism/StateValuesMTBDD.java
@@ -413,6 +413,30 @@ public class StateValuesMTBDD implements StateValues
 		return d;
 	}
 
+	@Override
+	public double maxFiniteOverBDD(JDDNode filter)
+	{
+		JDDNode tmp;
+		double d;
+
+		// filter filter
+		JDD.Ref(filter);
+		JDD.Ref(reach);
+		tmp = JDD.And(filter, reach);
+
+		// max of an empty set is -infinity
+		if (tmp.equals(JDD.ZERO)) return Double.NEGATIVE_INFINITY;
+
+		// set non-reach states to infinity
+		JDD.Ref(values);
+		tmp = JDD.ITE(tmp, values, JDD.MinusInfinity());
+
+		d = JDD.FindMaxFinite(tmp);
+		JDD.Deref(tmp);
+
+		return d;
+	}
+
 	@Override
 	public double sumOverBDD(JDDNode filter)
 	{
diff --git a/prism/src/prism/StateValuesVoid.java b/prism/src/prism/StateValuesVoid.java
index 5d7d81db..e231cdf7 100644
--- a/prism/src/prism/StateValuesVoid.java
+++ b/prism/src/prism/StateValuesVoid.java
@@ -178,6 +178,12 @@ public class StateValuesVoid implements StateValues
 		throw new UnsupportedOperationException();
 	}
 
+	@Override
+	public double maxFiniteOverBDD(JDDNode filter)
+	{
+		throw new UnsupportedOperationException();
+	}
+
 	@Override
 	public double sumOverBDD(JDDNode filter)
 	{