control.matlab.markov does not work

Hi,
I am new to this and it might be a silly doubt but since i could not find an answer anywhere, I am opening an issue here.

My problem is that I need to estimate initial states for my state space model. This requires that I find the markov parameters. I checked the documentation and tried the following example.
import control.matlab as mt inp=np.array([1,2]) outp=np.array([2,4]) mrk= mt.markov(outp,inp,1)

I just took one input and one output as shown above and tried to run the code. It always gives me the following error

"Traceback (most recent call last):

File "", line 3, in
mrk= mt.markov(outp,inp,1)

File "C:\Users\Gopal123\Anaconda3\lib\site-packages\control\modelsimp.py", line 421, in markov
Ulast = np.vstack((0, np.reshape(UU[0:n-1, m-2], (-1, 1))))

IndexError: too many indices for array
"
I am not able to understand the reason for this. I'd appreciate it if anyone can help me out. Thanks.

Here's the code from the unit tests for that function:

U = np.array([[1.], [1.], [1.], [1.], [1.]])
Y = U
M = 3
H = markov(Y,U,M)

This seems to work, though with a warning message.

@murrayrm

I think you found a bug in the markov function. Contrary to python-control convention, this function expects the time series of input and output data as a 2D array, where the time samples are by rows not columns.

Python 3.8.2 (default, Mar 05 2020, 18:58:42) [GCC]
Type "copyright", "credits" or "license" for more information.

IPython 7.13.0 -- An enhanced Interactive Python.

In [1]: import numpy as np

In [2]: from control.modelsimp import markov

In [3]: inp = np.array([[1,2]])

In [4]: outp = np.array([[2,4]])

In [5]: inp.T
Out[5]: 
array([[1],
       [2]])

In [6]: outp.T
Out[6]: 
array([[2],
       [4]])

In [7]: markov(outp,inp,1)
Traceback (most recent call last):

  File "<ipython-input-7-a693851d44f5>", line 1, in <module>
    markov(outp,inp,1)

  File "/usr/lib/python3.8/site-packages/control/modelsimp.py", line 424, in markov
    UU = np.hstack((UU, Ulast))

  File "<__array_function__ internals>", line 5, in hstack

  File "/usr/lib64/python3.8/site-packages/numpy/core/shape_base.py", line 344, in hstack
    return _nx.concatenate(arrs, 1)

  File "<__array_function__ internals>", line 5, in concatenate

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 1 and the array at index 1 has size 2


In [8]: markov(outp[:,0],inp[:,0],1)
Traceback (most recent call last):

  File "<ipython-input-8-319720f18954>", line 1, in <module>
    markov(outp[:,0],inp[:,0],1)

  File "/usr/lib/python3.8/site-packages/control/modelsimp.py", line 421, in markov
    Ulast = np.vstack((0, np.reshape(UU[0:n-1, m-2], (-1, 1))))

IndexError: too many indices for array


In [9]: markov(outp.T,inp.T,1)
/usr/lib/python3.8/site-packages/control/modelsimp.py:427: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  H = np.linalg.lstsq(UU, Y)[0]
Out[9]: 
array([[2.],
       [0.]])

In [10]: markov(outp.T,inp.T,2)
/usr/lib/python3.8/site-packages/control/modelsimp.py:427: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  H = np.linalg.lstsq(UU, Y)[0]
Out[10]: 
array([[2.],
       [0.]])

@murrayrm, if I fix this, would it break a lot of users' code? The function itself is 9 years old (d1b6b80) and the np.matrix to np.array change (de9e3e7) exposed the problem with 1D arrays.

Also, m<=2 does not make much sense. For state-space models, H[0], is just D.

@eigensteve

This may be a long shot, but I'd figure I try:

@eigensteve, the initial commit d1b6b80 looks like it might have been contributed by you 9 years ago. Could you please elaborate or point to some literature which algorithm is implemented here and what kind of input and output data it expects?

python-control/src/modelsimp.py

Lines 100 to 137 in d1b6b80

    
           def markov(Y,U,M): 
        
               """Calculate the first M Markov parameters [D CB CAB ...] from input U, output Y 
        
               Usage 
        
               ===== 
        
               H = markov(Y,U,M) 
        
               Currently only works for SISO 
        
               Inputs 
        
               ------ 
        
               Y : output data  
        
               U : input data 
        
               M : number of Markov parameters to output 
        
               Outputs 
        
               ------- 
        
               H : first M Markov parameters 
        
               """ 
        
               # Convert input parameters to matrices (if they aren't already) 
        
               Ymat = np.mat(Y) 
        
               Umat = np.mat(U) 
        
               n = np.size(U) 
        
               # Construct a matrix of control inputs to invert 
        
               UU = Umat 
        
               for i in range(1, M-1): 
        
                   newCol = np.vstack((0, UU[0:n-1,i-2])) 
        
                   UU = np.hstack((UU, newCol)) 
        
               Ulast = np.vstack((0, UU[0:n-1,M-2])) 
        
               for i in range(n-1,0,-1): 
        
                   Ulast[i] = np.sum(Ulast[0:i-1]) 
        
               UU = np.hstack((UU, Ulast)) 
        
               # Invert and solve for Markov parameters 
        
               H = UU.I 
        
               H = np.dot(H, Y)

The indices i-1 and M-2 in line 128 and 130 already raised some questions in de9e3e7.

My experiments with a randomly excited single mass oscillator did not result in any meaningful markov parameters.

Maybe a more modern and robust algorithm that can also deal with MIMO systems could replace the code?

Hi, Thank you for your response. Yes, it would be great if it can be used for MIMO systems as that is what I need it for. Also, it would be helpful if the documentation could be updated with the working example from above:
https://python-control.readthedocs.io/en/latest/generated/control.matlab.markov.html

Thanks for the email — it has been a while, but I do believe there are some better MIMO versions of the code. I will reach out to my lab and see if anyone has anything they can push. Best, Steve

…

On Apr 23, 2020, at 7:13 AM, Ben ***@***.***> wrote: This may be a long shot, but I'd figure I try: @eigensteve, the initial commit d1b6b80 looks like it might have been contributed by you 9 years ago. Could you please elaborate or point to some literature which algorithm is implemented here and what kind of input and output data it expects? https://github.com/python-control/python-control/blob/d1b6b809393e40748e71e9436d95d3313f11d104/src/modelsimp.py#L100-L137 The indices i-1 and M-2 in line 128 and 130 already raised some questions in de9e3e7. My experiments with a randomly excited single mass oscillator did not result in any meaningful markov parameters. Maybe a more modern and robust algorithm that can also deal with MIMO systems could replace the code? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

For release 0.8.4, I think we should leave the inconsistency in the time series convention (for backward compatibility) but perhaps check for this case and issue a warning message. One way to do this would be to add a keyword option that lets you switch around the convention and set the default to match the current use case.

For release 0.9.0, we can switch the default to match the time series convention but add to the code in the use_legacy_defaults (see PR #424) so that people can revert to the original convention if they have old code.

fixes to markov() + add tranpose keyword, default warning (#478)

Navigation Menu

Repository navigation

control.matlab.markov does not work #395

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	def markov(Y,U,M):
	"""Calculate the first M Markov parameters [D CB CAB ...] from input U, output Y

	Usage
	=====
	H = markov(Y,U,M)
	Currently only works for SISO

	Inputs
	------
	Y : output data
	U : input data
	M : number of Markov parameters to output

	Outputs
	-------
	H : first M Markov parameters

	"""

	# Convert input parameters to matrices (if they aren't already)
	Ymat = np.mat(Y)
	Umat = np.mat(U)
	n = np.size(U)

	# Construct a matrix of control inputs to invert
	UU = Umat
	for i in range(1, M-1):
	newCol = np.vstack((0, UU[0:n-1,i-2]))
	UU = np.hstack((UU, newCol))
	Ulast = np.vstack((0, UU[0:n-1,M-2]))
	for i in range(n-1,0,-1):
	Ulast[i] = np.sum(Ulast[0:i-1])
	UU = np.hstack((UU, Ulast))

	# Invert and solve for Markov parameters
	H = UU.I
	H = np.dot(H, Y)

control.matlab.markov does not work #395

Description

Activity

murrayrm commented on Apr 22, 2020

bnavigator commented on Apr 22, 2020

bnavigator commented on Apr 23, 2020

nrgopalrao commented on Apr 29, 2020

eigensteve commented on May 7, 2020

murrayrm commented on Jul 11, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions