au/transfer_function/
matrix.rs

1//! # Matrix of transfer functions
2//!
3//! `TfMatrix` allow the definition of a matrix of transfer functions. The
4//! numerators are stored in a matrix, while the denominator is stored once,
5//! since it is equal for every transfer function.
6//! * evaluation of a vector of inputs
7//! * conversion from a generic state space representation
8
9use ndarray::{Array2, Axis, Zip};
10use num_complex::Complex;
11use num_traits::{Float, MulAdd, One, Signed, Zero};
12
13use std::ops::{Index, IndexMut};
14use std::{
15    fmt,
16    fmt::{Debug, Display},
17};
18
19use crate::{
20    complex,
21    enums::Time,
22    linear_system::{self, SsGen},
23    polynomial::Poly,
24    polynomial_matrix::{MatrixOfPoly, PolyMatrix},
25};
26
27/// Matrix of transfer functions
28#[derive(Debug)]
29pub struct TfMatrix<T> {
30    /// Polynomial matrix of the numerators
31    num: MatrixOfPoly<T>,
32    /// Common polynomial denominator
33    den: Poly<T>,
34}
35
36/// Implementation of transfer function matrix
37impl<T> TfMatrix<T> {
38    /// Create a new transfer function matrix
39    ///
40    /// # Arguments
41    ///
42    /// * `num` - Polynomial matrix
43    /// * `den` - Characteristic polynomial of the system
44    fn new(num: MatrixOfPoly<T>, den: Poly<T>) -> Self {
45        Self { num, den }
46    }
47}
48
49impl<T: Clone> TfMatrix<T> {
50    /// Retrive the characteristic polynomial of the system.
51    #[must_use]
52    pub fn den(&self) -> Poly<T> {
53        self.den.clone()
54    }
55}
56
57impl<T: Float + MulAdd<Output = T>> TfMatrix<T> {
58    /// Evaluate the matrix transfers function.
59    ///
60    /// # Arguments
61    ///
62    /// * `s` - Array at which the Matrix of transfer functions is evaluated.
63    #[must_use]
64    pub fn eval(&self, s: &[Complex<T>]) -> Vec<Complex<T>> {
65        //
66        // ┌  ┐ ┌┌         ┐ ┌     ┐┐┌  ┐
67        // │y1│=││1/pc 1/pc│*│n1 n2│││s1│
68        // │y2│ ││1/pc 1/pc│ │n3 n4│││s2│
69        // └  ┘ └└         ┘ └     ┘┘└  ┘
70        // `*` is the element by element multiplication
71        // ┌     ┐ ┌┌         ┐ ┌     ┐┐ ┌┌     ┐ ┌     ┐┐
72        // │y1+y2│=││1/pc 1/pc│.│s1 s2││*││n1 n2│.│s1 s2││
73        // │y3+y4│ ││1/pc 1/pc│ │s1 s2││ ││n3 n4│ │s1 s2││
74        // └     ┘ └└         ┘ └     ┘┘ └└     ┘ └     ┘┘
75        // `.` means 'evaluated at'
76
77        // Create a matrix to contain the result of the evaluation.
78        let mut res = Array2::from_elem(
79            self.num.matrix().dim(),
80            Complex::<T>::new(T::zero(), T::zero()),
81        );
82
83        // Zip the result and the numerator matrix row by row.
84        Zip::from(res.genrows_mut())
85            .and(self.num.matrix().genrows())
86            .apply(|mut res_row, matrix_row| {
87                // Zip the row of the result matrix.
88                Zip::from(&mut res_row)
89                    .and(s) // The vector of the input.
90                    .and(matrix_row) // The row of the numerator matrix.
91                    .apply(|r, s, n| *r = complex::compdiv(n.eval(s), self.den.eval(s)));
92            });
93
94        res.sum_axis(Axis(1)).to_vec()
95    }
96}
97
98impl<T: Time> From<SsGen<f64, T>> for TfMatrix<f64> {
99    /// Convert a state-space representation into a matrix of transfer functions
100    ///
101    /// # Arguments
102    ///
103    /// `ss` - state space linear system
104    fn from(ss: SsGen<f64, T>) -> Self {
105        let (pc, a_inv) = linear_system::leverrier_f64(&ss.a);
106        let g = a_inv.left_mul(&ss.c).right_mul(&ss.b);
107        let rest = PolyMatrix::multiply(&pc, &ss.d);
108        let tf = g + rest;
109        Self::new(MatrixOfPoly::from(tf), pc)
110    }
111}
112
113impl<T: Time> From<SsGen<f32, T>> for TfMatrix<f32> {
114    /// Convert a state-space representation into a matrix of transfer functions
115    ///
116    /// # Arguments
117    ///
118    /// `ss` - state space linear system
119    fn from(ss: SsGen<f32, T>) -> Self {
120        let (pc, a_inv) = linear_system::leverrier_f32(&ss.a);
121        let g = a_inv.left_mul(&ss.c).right_mul(&ss.b);
122        let rest = PolyMatrix::multiply(&pc, &ss.d);
123        let tf = g + rest;
124        Self::new(MatrixOfPoly::from(tf), pc)
125    }
126}
127
128/// Implement read only indexing of the numerator of a transfer function matrix.
129///
130/// # Panics
131///
132/// Panics for out of bounds access.
133impl<T> Index<[usize; 2]> for TfMatrix<T> {
134    type Output = Poly<T>;
135
136    fn index(&self, i: [usize; 2]) -> &Poly<T> {
137        &self.num[i]
138    }
139}
140
141/// Implement mutable indexing of the numerator of a transfer function matrix
142/// returning its coefficients.
143///
144/// # Panics
145///
146/// Panics for out of bounds access.
147impl<T> IndexMut<[usize; 2]> for TfMatrix<T> {
148    fn index_mut(&mut self, i: [usize; 2]) -> &mut Poly<T> {
149        &mut self.num[i]
150    }
151}
152
153/// Implementation of transfer function matrix printing
154impl<T> fmt::Display for TfMatrix<T>
155where
156    T: Display + One + PartialEq + PartialOrd + Signed + Zero,
157{
158    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
159        let s_num = self.num.to_string();
160        let s_den = self.den.to_string();
161
162        let length = s_den.len();
163        let dash = "\u{2500}".repeat(length);
164
165        write!(f, "{}\n{}\n{}", s_num, dash, s_den)
166    }
167}
168
169#[cfg(test)]
170mod tests {
171    use super::*;
172    use crate::{poly, Ss, Ssd};
173
174    #[test]
175    fn tf_matrix_new() {
176        let sys = Ssd::new_from_slice(
177            2,
178            2,
179            2,
180            &[-2., 0., 0., -1.],
181            &[0., 1., 1., 2.],
182            &[1., 2., 3., 4.],
183            &[1., 0., 0., 1.],
184        );
185        let tfm = TfMatrix::<f32>::from(sys);
186        assert_eq!(tfm[[0, 0]], poly!(6., 5., 1.));
187        assert_eq!(tfm[[0, 1]], poly!(9., 5.));
188        assert_eq!(tfm[[1, 0]], poly!(8., 4.));
189        assert_eq!(tfm[[1, 1]], poly!(21., 14., 1.));
190        assert_eq!(tfm.den(), poly!(2., 3., 1.));
191    }
192
193    #[test]
194    fn tf_matrix_eval() {
195        let sys = Ss::new_from_slice(
196            2,
197            2,
198            2,
199            &[-2., 0., 0., -1.],
200            &[0., 1., 1., 2.],
201            &[1., 2., 3., 4.],
202            &[1., 0., 0., 1.],
203        );
204        let tfm = TfMatrix::from(sys);
205        let i = Complex::<f64>::i();
206        let res = tfm.eval(&[i, i]);
207        assert_relative_eq!(res[0].re, 4.4, max_relative = 1e-15);
208        assert_relative_eq!(res[0].im, -3.2, max_relative = 1e-15);
209        assert_relative_eq!(res[1].re, 8.2, max_relative = 1e-15);
210        assert_relative_eq!(res[1].im, -6.6, max_relative = 1e-15);
211    }
212
213    #[test]
214    fn tf_matrix_index_mut() {
215        let sys = Ss::new_from_slice(
216            2,
217            2,
218            2,
219            &[-2., 0., 0., -1.],
220            &[0., 1., 1., 2.],
221            &[1., 2., 3., 4.],
222            &[1., 0., 0., 1.],
223        );
224        let mut tfm = TfMatrix::from(sys);
225        tfm[[0, 0]] = poly!(1., 2., 3.);
226        assert_eq!(poly!(1., 2., 3.), tfm[[0, 0]]);
227    }
228
229    #[test]
230    fn tf_matrix_print() {
231        let sys = Ss::new_from_slice(
232            2,
233            2,
234            2,
235            &[-2., 0., 0., -1.],
236            &[0., 1., 1., 2.],
237            &[1., 2., 3., 4.],
238            &[1., 0., 0., 1.],
239        );
240        let tfm = TfMatrix::from(sys);
241        assert_eq!(
242            "[[6 +5s +1s^2, 9 +5s],\n [8 +4s, 21 +14s +1s^2]]\n\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\n2 +3s +1s^2",
243            format!("{}", tfm)
244        );
245    }
246}
au/transfer_function/matrix.rs

au/transfer_function/
matrix.rs
