g-bauer · Oct 31, 2026 · Oct 29, 2026
diff --git a/crates/feos-core/src/phase_equilibria/stability_analysis.rs b/crates/feos-core/src/phase_equilibria/stability_analysis.rs
@@ -16,6 +16,7 @@ const MIN_EIGENVAL: f64 = 1E-03;
 const ETA_STEP: f64 = 0.25;
 const MINIMIZE_KMAX: usize = 100;
 const ZERO_TPD: f64 = -1E-08;
+const DUPLICATE_TOL: f64 = 1E-03;
 
 /// # Stability analysis
 impl<E: Residual<N>, N: Gradients> State<E, N>
@@ -31,37 +32,106 @@ where
     /// Perform a stability analysis. The result is a list of [State]s with
     /// negative tangent plane distance (i.e. lower Gibbs energy) that can be
     /// used as initial estimates for a phase equilibrium calculation.
+    /// The states are sorted by their tangent plane distance, starting with
+    /// the most negative one.
+    ///
+    /// A trial phase whose minimization fails is still returned if its
+    /// tangent plane distance is negative, because that already proves the
+    /// instability. Otherwise it is skipped. An error is only returned if
+    /// no trial phase could be evaluated successfully.
     pub fn stability_analysis(&self, options: SolverOptions) -> FeosResult<Vec<State<E, N>>> {
-        let mut result = Vec::new();
+        let mut result: Vec<(f64, State<E, N>)> = Vec::new();
+        let mut successful_trials = 0;
+        let mut last_error = None;
         for i_trial in 0..self.eos.components() + 1 {
             let phase = if i_trial == self.eos.components() {
                 "Vapor phase".to_string()
             } else {
                 format!("Liquid phase {}", i_trial + 1)
             };
-            if let Ok(mut trial_state) = self.define_trial_state(i_trial) {
-                let (tpd, i) = self.minimize_tpd(&mut trial_state, options)?;
-                let msg = if let Some(tpd) = tpd {
-                    if tpd < ZERO_TPD {
-                        if result
-                            .iter()
-                            .any(|s| PhaseEquilibrium::is_trivial_solution(s, &trial_state))
-                        {
-                            "Found already identified minimum"
-                        } else {
-                            result.push(trial_state);
-                            "Found candidate"
-                        }
+            let mut trial_state = match self.define_trial_state(i_trial) {
+                Ok(trial_state) => trial_state,
+                Err(e) => {
+                    log_result!(options.verbosity, "{}: Failed ({})\n", phase, e);
+                    last_error = Some(e);
+                    continue;
+                }
+            };
+            // The tpd tracked during the minimization is only an estimate, so the
+            // true tpd of the final trial state is used to compare candidates.
+            let (tpd, steps) = match self.minimize_tpd(&mut trial_state, options) {
+                Ok((None, i)) => (None, Some(i)),
+                Ok((Some(_), i)) => (Some(self.tangent_plane_distance(&trial_state)), Some(i)),
+                Err(e) => {
+                    let tpd = self.tangent_plane_distance(&trial_state);
+                    if tpd < ZERO_TPD && !PhaseEquilibrium::is_trivial_solution(self, &trial_state)
+                    {
+                        (Some(tpd), None)
+                    } else {
+                        log_result!(options.verbosity, "{}: Failed ({})\n", phase, e);
+                        last_error = Some(e);
+                        continue;
+                    }
+                }
+            };
+            successful_trials += 1;
+            let msg = if let Some(tpd) = tpd {
+                if tpd < ZERO_TPD {
+                    if result.iter().any(|(_, s)| s.is_same_phase(&trial_state)) {
+                        "Found already identified minimum"
                     } else {
-                        "Found minimum > 0"
+                        result.push((tpd, trial_state));
+                        "Found candidate"
                     }
                 } else {
-                    "Found trivial solution"
-                };
-                log_result!(options.verbosity, "{}: {} in {} step(s)\n", phase, msg, i);
+                    "Found minimum > 0"
+                }
+            } else {
+                "Found trivial solution"
+            };
+            match steps {
+                Some(i) => {
+                    log_result!(options.verbosity, "{}: {} in {} step(s)\n", phase, msg, i)
+                }
+                None => log_result!(options.verbosity, "{}: {} (not converged)\n", phase, msg),
             }
         }
-        Ok(result)
+        if let Some(e) = last_error
+            && successful_trials == 0
+        {
+            return Err(e);
+        }
+        result.sort_by(|(tpd1, _), (tpd2, _)| tpd1.total_cmp(tpd2));
+        Ok(result.into_iter().map(|(_, s)| s).collect())
+    }
+
+    /// Check if two candidates from the stability analysis describe the same phase.
+    ///
+    /// In contrast to [PhaseEquilibrium::is_trivial_solution], the partial densities
+    /// are compared relative to the total density, so that differences in trace
+    /// components (which can be large in relative terms) are ignored.
+    fn is_same_phase(&self, other: &State<E, N>) -> bool {
+        let rho1 = self.density.into_reduced();
+        let rho2 = other.density.into_reduced();
+        let max_diff = self
+            .molefracs
+            .zip_fold(&other.molefracs, 0.0, |acc: f64, x1, x2| {
+                acc.max((x1 * rho1 - x2 * rho2).abs())
+            });
+        max_diff / rho1.max(rho2) < DUPLICATE_TOL
+    }
+
+    /// Tangent plane distance of the trial state with respect to this state.
+    ///
+    /// Components that are absent in the trial state do not contribute (w ln w -> 0).
+    pub(super) fn tangent_plane_distance(&self, trial_state: &State<E, N>) -> f64 {
+        let ln_phi_z = self.ln_phi();
+        let ln_phi_w = trial_state.ln_phi();
+        let (z, w) = (&self.molefracs, &trial_state.molefracs);
+        (0..w.len())
+            .filter(|&i| w[i] > 0.0)
+            .map(|i| w[i] * (w[i].ln() + ln_phi_w[i] - z[i].ln() - ln_phi_z[i]))
+            .sum()
     }
 
     fn define_trial_state(&self, dominant_component: usize) -> FeosResult<State<E, N>> {

diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs
@@ -248,14 +248,6 @@ where
 
         Ok(new_vle_state)
     }
-
-    fn tangent_plane_distance(&self, trial_state: &State<E, N>) -> f64 {
-        let ln_phi_z = self.ln_phi();
-        let ln_phi_w = trial_state.ln_phi();
-        let z = &self.molefracs;
-        let w = &trial_state.molefracs;
-        w.dot(&(w.map(f64::ln) + ln_phi_w - z.map(f64::ln) - ln_phi_z))
-    }
 }
 
 impl<E: Residual<N>, N: Gradients> PhaseEquilibrium<E, 2, N>
@@ -420,9 +412,12 @@ where
     }
 
     fn vle_init_stability(feed_state: &State<E, N>) -> FeosResult<(Self, Option<Self>)> {
-        let mut stable_states = feed_state.stability_analysis(SolverOptions::default())?;
-        let state1 = stable_states.pop();
-        let state2 = stable_states.pop();
+        // candidates are sorted by tpd, so the most negative ones are used
+        let mut stable_states = feed_state
+            .stability_analysis(SolverOptions::default())?
+            .into_iter();
+        let state1 = stable_states.next();
+        let state2 = stable_states.next();
         if let Some(s1) = state1 {
             let init1 = if s1.density < feed_state.density {
                 Self::two_phase(s1.clone(), feed_state.clone())

diff --git a/crates/feos/tests/pcsaft/stability_analysis.rs b/crates/feos/tests/pcsaft/stability_analysis.rs
@@ -1,6 +1,7 @@
 use feos::pcsaft::{PcSaft, PcSaftParameters};
 use feos_core::parameter::IdentifierOption;
-use feos_core::{DensityInitialization, PhaseEquilibrium, SolverOptions, State};
+use feos_core::{DensityInitialization, PhaseEquilibrium, Residual, SolverOptions, State};
+use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
 
@@ -48,3 +49,168 @@ fn test_stability_analysis() -> Result<(), Box<dyn Error>> {
     assert!(liquid_check.is_empty());
     Ok(())
 }
+
+#[test]
+fn test_stability_analysis_failed_trial() -> Result<(), Box<dyn Error>> {
+    let params = PcSaftParameters::from_json(
+        vec!["water_np", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?;
+    let mix = PcSaft::new(params);
+    let unstable = State::new_npt(
+        &&mix,
+        300.0 * KELVIN,
+        1.0 * BAR,
+        0.5,
+        Some(DensityInitialization::Liquid),
+    )?;
+
+    // the second liquid trial phase needs more than 4 iterations; the
+    // remaining trial phases still find the instability
+    let options = SolverOptions {
+        max_iter: Some(4),
+        ..Default::default()
+    };
+    let check = unstable.stability_analysis(options)?;
+    assert!(!check.is_empty());
+
+    // after a single iteration, the minimizations are not converged, but the
+    // negative tangent plane distances already prove the instability
+    let options = SolverOptions {
+        max_iter: Some(1),
+        ..Default::default()
+    };
+    let check = unstable.stability_analysis(options)?;
+    assert_eq!(check.len(), 2);
+
+    // if the minimization fails for all trial phases without finding a
+    // negative tangent plane distance, the error is returned
+    let params = PcSaftParameters::from_json(
+        vec!["propane", "butane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?;
+    let mix = PcSaft::new(params);
+    let stable = State::new_npt(
+        &&mix,
+        300.0 * KELVIN,
+        20.0 * BAR,
+        0.5,
+        Some(DensityInitialization::Liquid),
+    )?;
+    assert!(stable.stability_analysis(options).is_err());
+    Ok(())
+}
+
+#[test]
+fn test_stability_analysis_sorted() -> Result<(), Box<dyn Error>> {
+    let params = PcSaftParameters::from_json(
+        vec!["water_np", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?;
+    let mix = PcSaft::new(params);
+    let unstable = State::new_npt(
+        &&mix,
+        350.0 * KELVIN,
+        1.0 * BAR,
+        0.5,
+        Some(DensityInitialization::Liquid),
+    )?;
+
+    // the trial phases find a water-rich liquid, a hexane-rich liquid and a
+    // vapor (in that order); the vapor has the most negative tpd
+    let check = unstable.stability_analysis(SolverOptions::default())?;
+    assert_eq!(check.len(), 3);
+    let tpds: Vec<_> = check.iter().map(|s| tpd(&unstable, s)).collect();
+    assert!(tpds.is_sorted());
+    assert!(check[0].density < check[1].density);
+    assert!(check[0].density < check[2].density);
+
+    // converged and not converged candidates are sorted by the same measure
+    let unstable = State::new_npt(
+        &&mix,
+        300.0 * KELVIN,
+        5.0 * BAR,
+        0.1,
+        Some(DensityInitialization::Vapor),
+    )?;
+    let options = SolverOptions {
+        max_iter: Some(3),
+        ..Default::default()
+    };
+    let check = unstable.stability_analysis(options)?;
+    assert_eq!(check.len(), 2);
+    let tpds: Vec<_> = check.iter().map(|s| tpd(&unstable, s)).collect();
+    assert!(tpds.is_sorted());
+    Ok(())
+}
+
+#[test]
+fn test_stability_analysis_zero_component() -> Result<(), Box<dyn Error>> {
+    let params = PcSaftParameters::from_json(
+        vec!["propane", "butane", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?;
+    let mix = PcSaft::new(params);
+    let unstable = State::new_npt(
+        &&mix,
+        300.0 * KELVIN,
+        1.2 * BAR,
+        &dvector![0.0, 0.5, 0.5],
+        Some(DensityInitialization::Liquid),
+    )?;
+
+    // absent components must not turn the tpd into NaN, also for
+    // trial phases that are not converged
+    for max_iter in [1, 100] {
+        let options = SolverOptions {
+            max_iter: Some(max_iter),
+            ..Default::default()
+        };
+        let check = unstable.stability_analysis(options)?;
+        assert!(!check.is_empty());
+        assert!(check.iter().all(|s| tpd(&unstable, s) < 0.0));
+    }
+    Ok(())
+}
+
+#[test]
+fn test_stability_analysis_duplicates() -> Result<(), Box<dyn Error>> {
+    let params = PcSaftParameters::from_json(
+        vec!["water_np", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?;
+    let mix = PcSaft::new(params);
+    let unstable = State::new_npt(
+        &&mix,
+        300.0 * KELVIN,
+        20.0 * BAR,
+        0.9,
+        Some(DensityInitialization::Liquid),
+    )?;
+
+    // the hexane-rich liquid is found from two trial phases; it only differs
+    // in the (trace) water content and must be reported once
+    let check = unstable.stability_analysis(SolverOptions::default())?;
+    assert_eq!(check.len(), 2);
+    Ok(())
+}
+
+fn tpd<E: Residual>(feed: &State<E>, trial: &State<E>) -> f64 {
+    let d = feed.molefracs.map(f64::ln) + feed.ln_phi();
+    let ln_phi = trial.ln_phi();
+    let w = &trial.molefracs;
+    (0..w.len())
+        .filter(|&i| w[i] > 0.0)
+        .map(|i| w[i] * (w[i].ln() + ln_phi[i] - d[i]))
+        .sum()
+}