Bugs/norms dim zero

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "control_systems_torbox"
-version = "0.2.0"
+version = "0.2.1"
 edition = "2024"
 authors = ["Tor Børve Rasmussen <tor.jobb.skule@gmail.com>"]
 description = "Control systems toolbox"

build.rs

@@ -13,7 +13,7 @@ fn main() {
         .to_string();
     let make_args: Vec<String> = vec![
         "-j8".to_string(),
-        format!("BUILD_DIR={}", slicot_build_dir).to_string(),
+        format!("BUILD_DIR={slicot_build_dir}").to_string(),
     ];
 
     let status = make()
@@ -32,7 +32,7 @@ fn main() {
         panic!("Failed to build SLICOT")
     }
 
-    println!("cargo:rustc-link-search={}", slicot_build_dir);
+    println!("cargo:rustc-link-search={slicot_build_dir}");
     println!("cargo:rustc-link-lib=static=slicot");
     println!("cargo:rustc-link-lib=static=lpkaux");
     println!("cargo:rustc-link-lib=gfortran");

src/analysis/frequency_response.rs

@@ -254,8 +254,7 @@ fn freq_response_ss_mat(
 
     if info != 0 {
         return Err(format!(
-            "Failed to compute frequency response of state space system. Slicot TB05AD failed with error code: {}",
-            info
+            "Failed to compute frequency response of state space system. Slicot TB05AD failed with error code: {info}"
         ));
     }
 

src/analysis/system_properties.rs

@@ -50,6 +50,17 @@ fn h2_norm<U: Time + 'static>(
     let m = b.ncols();
     let p = c.nrows();
 
+    assert!(m > 0);
+    assert!(p > 0);
+
+    if n == 0 {
+        let d_all_zero = d.iter().all(|&x| x == 0.0);
+        if d_all_zero {
+            return Ok(0.);
+        }
+        return Err("Order of system is zero and D matrix is not zero".into());
+    }
+
     let lda = a.nrows();
     let ldb = b.nrows();
     let ldc = c.nrows();
@@ -66,6 +77,7 @@ fn h2_norm<U: Time + 'static>(
     let mut iwarn = -1 as c_int;
     let mut info = -1 as c_int;
 
+    // println!("lda: {}", lda);
     let h2_norm = unsafe {
         ab13bd_(
             time_domain.as_ptr(),
@@ -92,13 +104,11 @@ fn h2_norm<U: Time + 'static>(
 
     if info != 0 {
         return Err(Box::new(std::io::Error::other(format!(
-            "SLICOT ab13bd failed with the follwing error indicator: {}",
-            info
+            "SLICOT ab13bd failed with the follwing error indicator: {info}"
         ))));
     } else if iwarn != 0 {
         return Err(Box::new(std::io::Error::other(format!(
-            "SLICOT ab13bd returned with warning about numerical stability. Iwarn = {}",
-            iwarn
+            "SLICOT ab13bd returned with warning about numerical stability. Iwarn = {iwarn}"
         ))));
     }
 
@@ -122,6 +132,21 @@ fn hinf_norm<U: Time + 'static>(
 ) -> Result<f64, String> {
     Ss::<U>::verify_dimensions(&a, &b, &c, &d).map_err(|e| e.to_string())?;
 
+    let n = a.nrows();
+    let m = b.ncols();
+    let p = c.nrows();
+
+    assert!(p > 0);
+    assert!(m > 0);
+
+    if n == 0 {
+        // The gain is constant and equal to D
+        // H-inf norm is the largest singular value
+        let singular_values = d.singular_values();
+        let max_singular_value = singular_values.max();
+        return Ok(max_singular_value);
+    }
+
     let time_domain = if std::any::TypeId::of::<U>()
         == std::any::TypeId::of::<Continuous>()
     {
@@ -137,10 +162,6 @@ fn hinf_norm<U: Time + 'static>(
     let equilibration = CString::new("S").unwrap(); // Perform scaling
     let d_is_nonzero = CString::new("D").unwrap();
 
-    let n = a.nrows();
-    let m = b.ncols();
-    let p = c.nrows();
-
     let mut e = DMatrix::identity(n, n);
 
     let mut f_peak = [0.0, 1.0]; // initial guess not active
@@ -199,8 +220,7 @@ fn hinf_norm<U: Time + 'static>(
 
     if info != 0 {
         return Err(format!(
-            "SLICOT ab13dd returned with error code. info = {}",
-            info
+            "SLICOT ab13dd returned with error code. info = {info}"
         ));
     }
 
@@ -293,8 +313,7 @@ pub fn zeros(
 
     if info != 0 {
         return Err(format!(
-            "SLICOT failed to find zeros of state-space system. Info = {}",
-            info
+            "SLICOT failed to find zeros of state-space system. Info = {info}"
         ));
     }
     let ldwork = dwork[0] as usize;
@@ -336,8 +355,7 @@ pub fn zeros(
 
     if info != 0 {
         return Err(format!(
-            "SLICOT failed to find zeros of state-space system. Info = {}",
-            info
+            "SLICOT failed to find zeros of state-space system. Info = {info}"
         ));
     }
     if num_inv_zeros == 0 {
@@ -384,7 +402,7 @@ pub fn zeros(
         );
     }
     if info != 0 {
-        return Err(format!("LAPACK DGGEV returned info = {}", info));
+        return Err(format!("LAPACK DGGEV returned info = {info}"));
     }
     lwork = dwork[0] as c_int;
     assert!(lwork > 0);
@@ -413,7 +431,7 @@ pub fn zeros(
         );
     }
     if info != 0 {
-        return Err(format!("LAPACK DGGEV returned info = {}", info));
+        return Err(format!("LAPACK DGGEV returned info = {info}"));
     }
 
     let mut zeros = Vec::with_capacity(n_gen_eigen);
@@ -446,6 +464,7 @@ impl<U: Time + 'static> Ss<U> {
     /// - `Ok(f64)`: The H2 norm of the system.
     /// - `Err(String)`: An error message if the computation fails.
     pub fn norm_h2(&self) -> Result<f64, String> {
+        // println!("before h2: {}", self);
         h2_norm::<U>(
             self.a().clone(),
             self.b().clone(),
@@ -510,6 +529,7 @@ impl<U: Time + 'static> Ss<U> {
 #[cfg(test)]
 mod tests {
     use crate::{
+        Continuous, Ss,
         systems::Tf,
         transformations::SsRealization::{ControllableCF, ObservableCF},
     };
@@ -546,6 +566,13 @@ mod tests {
                 epsilon = 1e-2
             );
         }
+
+        let sys = Ss::<Continuous>::new_from_scalar(-2.0);
+        assert_abs_diff_eq!(sys.norm_hinf().unwrap(), 2.0);
+
+        assert!(sys.norm_h2().is_err());
+        let sys = Ss::<Continuous>::new_from_scalar(0.0);
+        assert_eq!(sys.norm_h2().unwrap(), 0.0);
     }
 
     #[test]
@@ -592,9 +619,9 @@ mod tests {
     #[test]
     fn ss_zeros() {
         let tf = (Tf::s() - 1.0) * (Tf::s() + 4.0) / (Tf::s() + 2.0).powi(2);
-        println!("tf: \n{}", tf);
+        println!("tf: \n{tf}");
         let zeros = tf.to_ss_method(ControllableCF).unwrap().zeros().unwrap();
-        println!("zeros: {:?}", zeros);
+        println!("zeros: {zeros:?}");
         assert_eq!(zeros.len(), 2);
 
         let tf = 1.0 / Tf::s();
@@ -643,18 +670,18 @@ mod tests {
     #[test]
     fn ss_is_stable() {
         let tf = 1.0 / Tf::s();
-        assert_eq!(tf.to_ss_method(ObservableCF).unwrap().is_stable(), false);
+        assert!(!tf.to_ss_method(ObservableCF).unwrap().is_stable());
 
         let tf = (Tf::s() + 1.0) / (Tf::s() + 1.0).powi(4);
         let ss = tf.to_ss_method(ObservableCF).unwrap();
-        assert_eq!(ss.is_stable(), true);
+        assert!(ss.is_stable());
 
         let tf = 1.0 / ((Tf::s() + 1.0) * (Tf::s() - 2.0));
         let ss = tf.to_ss_method(ObservableCF).unwrap();
-        assert_eq!(ss.is_stable(), false);
+        assert!(!ss.is_stable());
 
         let tf = 1.0 / (Tf::s().powi(2) + 2.0 * 0.01 * Tf::s() + 1.0);
         let ss = tf.to_ss_method(ObservableCF).unwrap();
-        assert_eq!(ss.is_stable(), true);
+        assert!(ss.is_stable());
     }
 }

src/systems/polynom.rs

@@ -60,8 +60,8 @@ where
 
             let term_i = match exp {
                 0 => val.to_string(),
-                1 => format!("{}{}", val, var_name),
-                _ => format!("{}{}^{}", val, var_name, exp),
+                1 => format!("{val}{var_name}"),
+                _ => format!("{val}{var_name}^{exp}"),
             };
             terms.push(term_i);
         }
@@ -76,7 +76,7 @@ where
 {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         let str_poly = self.to_string_variable("x");
-        write!(f, "{}", str_poly)
+        write!(f, "{str_poly}")
     }
 }
 
@@ -353,11 +353,7 @@ where
         assert!(num_str.len() <= max_len);
         assert!(den_str.len() <= max_len);
 
-        write!(
-            formatter,
-            "\n  {}\n  {}\n  {}\n\n",
-            num_str, dash_line, den_str
-        )
+        write!(formatter, "\n  {num_str}\n  {dash_line}\n  {den_str}\n\n")
     }
 }
 
@@ -649,7 +645,7 @@ mod tests {
             lhs_vals
                 .iter()
                 .zip(rhs_vals.iter())
-                .all(|(x, y)| f64::abs_diff_eq(&x, &y, epsilon))
+                .all(|(x, y)| f64::abs_diff_eq(x, y, epsilon))
         }
     }
 
@@ -712,9 +708,9 @@ mod tests {
         }
 
         {
-            println!("poly: {}", p);
+            println!("poly: {p}");
             let p_str = p.to_string();
-            println!("poly to string: {}", p_str);
+            println!("poly to string: {p_str}");
         }
     }
 
@@ -760,7 +756,7 @@ mod tests {
         assert_abs_diff_eq!(rf.num.coeffs[0], 1.0);
         assert_abs_diff_eq!(rf.num.coeffs.last().unwrap().clone(), 3.0);
 
-        println!("Rational function normalize: \n{}", rf);
+        println!("Rational function normalize: \n{rf}");
     }
 
     #[test]

src/systems/state_space.rs

@@ -612,7 +612,7 @@ impl<U: Time + 'static> fmt::Display for Ss<U> {
         } else {
             "Unknown"
         };
-        write!(f, "{}-time state-space model\n\n", time_domain)
+        write!(f, "{time_domain}-time state-space model\n\n")
     }
 }
 
@@ -635,7 +635,7 @@ mod tests {
         rng: &mut U,
         max_order: usize,
     ) -> Tf<f64, Continuous> {
-        let den_order = rng.random_range(1..=max_order) as usize;
+        let den_order = rng.random_range(1..=max_order);
         let num_order = rng.random_range(0..=den_order);
 
         let num: Vec<f64> = (0..=num_order)
@@ -884,14 +884,14 @@ mod tests {
     #[test]
     fn ss_display() {
         let sys = (1.0 / Tf::s()).to_ss().unwrap();
-        println!("1: {}", sys);
+        println!("1: {sys}");
 
         let sys = Ss::<Discrete>::new_from_scalar(2.0);
-        println!("2: {}", sys);
+        println!("2: {sys}");
 
         let sys = ((1.0 + Tf::s()) / ((Tf::s() - 1.0) * Tf::s()).powi(2))
             .to_ss()
             .unwrap();
-        println!("3: {}after", sys);
+        println!("3: {sys}after");
     }
 }

src/systems/transfer_function.rs

@@ -497,10 +497,10 @@ mod tests {
             assert_abs_diff_eq!(y.re, y_expect.re);
             assert_abs_diff_eq!(y.im, y_expect.im);
         }
-        println!("Tf cont: \n{}", tf);
+        println!("Tf cont: \n{tf}");
 
         let tf_z = Tf::<f64, Discrete>::new_from_rf(tf.rf);
-        println!("Tf discrete: \n {}", tf_z);
+        println!("Tf discrete: \n {tf_z}");
     }
 
     #[test]
@@ -553,16 +553,16 @@ mod tests {
             Tf::<f64, Discrete>::new(&[0.0, 0., 0., 1., 2., 0.0], &[0.0, 1.]);
         assert_eq!(sys.degree_num_den(), (4, 1));
         assert_eq!(sys.relative_degree(), 4 - 1);
-        assert_eq!(sys.is_proper(), false);
-        assert_eq!(sys.is_strictly_proper(), false);
+        assert!(!sys.is_proper());
+        assert!(!sys.is_strictly_proper());
 
         let sys = Tf::s() / Tf::s();
-        assert_eq!(sys.is_proper(), true);
-        assert_eq!(sys.is_strictly_proper(), false);
+        assert!(sys.is_proper());
+        assert!(!sys.is_strictly_proper());
 
         let sys = 1. / Tf::s();
-        assert_eq!(sys.is_proper(), true);
-        assert_eq!(sys.is_strictly_proper(), true);
+        assert!(sys.is_proper());
+        assert!(sys.is_strictly_proper());
     }
 
     #[test]
@@ -582,11 +582,11 @@ mod tests {
     #[test]
     fn tf_display() {
         let tf = 1.0 / Tf::s();
-        println!("1: {}", tf);
+        println!("1: {tf}");
         let tf = 1.0 / Tf::s().powi(5) * (Tf::s() + 1.0);
-        println!("2: {}", tf);
+        println!("2: {tf}");
         let tf =
             2.1 / Tf::s() * (Tf::s() - 1.2).powi(5) / (Tf::s() + 1.0).powi(3);
-        println!("3: {}", tf);
+        println!("3: {tf}");
     }
 }