Re: matrix exponential used by the control package

On Tue, 2020-07-21 at 07:22 +0200, Torsten Lilge wrote:
> On Mon, 2020-07-20 at 17:22 -0400, Doug Stewart wrote:
> > 
> > On Mon, Jul 20, 2020 at 4:38 PM Torsten Lilge1960 Blue Heron Drive,
> > <
> > [email protected]> wrote:
> > > Hi everyone,
> > > 
> > > the control package uses the SLICOT function mb05nd [1]
> > > (__sl_mb05nd__()) instead of the octave function expm() for the
> > > matrix
> > > exponential. The problem with mb05nd is that when discretizing a
> > > static
> > > gain in state space results in 
> > > 
> > > > > [a,b,c,d] = ssdata (c2d (ss (tf (1,1)),1))
> > > a = [](1x0)
> > > b = 0
> > > c = [](1x0)
> > > d = 1
> > > 
> > > Instead, using expm() leads to the correct result
> > > 
> > > > > [a,b,c,d] = ssdata (c2d (ss (tf (1,1)),1))
> > > a = [](0x0)
> > > b = [](0x1)
> > > c =
> > > [](1x0)
> > > d = 1
> > > 
> > > Are there any reasons against using expm()?
> > > 
> > > Torsten
> > > 
> > > 
> > > [1] http://slicot.org/objects/software/shared/doc/MB05ND.html
> > > 
> > > 
> > 
> > Is there a difference in speed?
> > If it is slower then we should look and see where the '0" is coming
> > from.
> 
> __sl_mb05nd__ is a compiled function and - at least on my system -
> approx. 20 times faster.
> 

Hello Doug,

The routine mb05nd does not accept input dimensions 0x0 for an empty
matrix but it accepts (1x0) which is forced by the wrapper sl_mb05nd for
empty matrices. I have therefore  pushed changeset 7a8204f62fe1 to the
control package which adds the following test to sl_mb05nd:

+        if (a.numel () == 0)
+          {
+            // given matrix is empty, just return the empty matrix,
+            // otherise mb05d would make a [](1x0) matrix from it
+            retval(0) = a;
+            retval(1) = a;
+            return retval;
+          }

I have also prepared a patch for adding the first order-hold method
'foh' to c2d. This enables us to use 'foh' in lsim() or ramp(). Should I
push this to the default branch or should I create a feature branch for
this?

Torsten