Issue with Non-Zero starting conditions when using control.forced_response - Not working as expected? #843

TheWicklowWolf · 2023-01-05T16:11:56Z

TheWicklowWolf
Jan 5, 2023

When using forced_response with zero initial conditions it works as expected:

#Sample Data
t=[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
 54, 55, 56, 57, 58, 59, 60]

pv=[13.45,13.45,13.57,13.64,13.43,13.63,13.5,13.63,13.5,14.11,15.13,16.04
,16.88,17.84,18.64,19.2,19.95,20.56,21.33,21.9,22.49,22.86,23.38,23.97
,24.47,24.79,25.14,25.62,26.15,26.28,26.62,26.94,27.52,27.62,28.11,28.37
,28.36,28.87,29.1,29.16,29.55,29.66,29.86,30.08,30.13,30.19,30.52,30.7
,30.82,30.99,31.13,31.2,31.29,31.31,31.5,31.5,31.61,31.82,31.79,31.84
,31.96]

cv=[ 0,  0,  0,  0,  0,  0,  0,  0,  0, 10, 10, 10, 10, 10, 10, 10, 10, 10,
 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
 10, 10, 10, 10, 10, 10, 10]

import numpy as np
import matplotlib.pyplot as plt
import control

gain=2
tc=20
bias=13.5

def fo_model(gain, tc):
    num=np.array([gain])
    den=np.array([tc,1])
    H = control.tf(num,den)
    (t_ ,y) = control.forced_response(H, t, cv, X0=0)
    return y+bias

plt.plot(t, pv, 'b', label='PV')
plt.plot(t, cv, 'r', label='CV')
plt.plot(t, fo_model(gain,tc), 'g', label='FO')
plt.legend()
plt.show()

CV goes from 0-10 and produces this, which is as expected:

However it doesn't work correctly (or rather it is not what I expected) when initial conditions are not zero:

#Slightly different Sample Data
t=[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
 54, 55, 56, 57, 58, 59, 60]

pv=[33.53,33.37,33.57,33.35,33.53,33.36,33.52,33.31,33.57,33.86,34.69,35.54
,36.61,37.33,38.08,38.75,39.61,40.37,40.91,41.57,42.01,42.6,43.28,43.78
,44.16,44.58,45.03,45.48,45.9,46.07,46.69,46.92,47.16,47.6,47.92,48.01
,48.26,48.49,48.94,48.97,49.39,49.39,49.67,49.8,50.02,50.15,50.5,50.42
,50.68,50.75,51.01,51.06,51.25,51.29,51.44,51.48,51.68,51.57,51.89,51.96
,51.91]

cv=[10, 10, 10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 20, 20, 20, 20,
 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
 20, 20, 20, 20, 20, 20, 20]

import numpy as np
import matplotlib.pyplot as plt
import control

gain=2
tc=20
bias=13.5

def fo_model(gain, tc):
    num=np.array([gain])
    den=np.array([tc,1])
    H = control.tf(num,den)
    (t_ ,y) = control.forced_response(H, t, cv, X0=pv[0])
    return y+bias

plt.plot(t, pv, 'b', label='PV')
plt.plot(t, cv, 'r', label='CV')
plt.plot(t, fo_model(gain,tc), 'g', label='FO')
plt.legend()
plt.show()

CV now goes from 10-20 (same system so it should match) but it produces:

Any ideas about where I'm going wrong?
Just found this library and it's excellent, I appreciate the work!

Answered by bnavigator

Jan 5, 2023

The code you posted should have given you a UserWarning:

UserWarning: Non-zero initial condition given for transfer function system. Internal conversion to state space used; may not be consistent with given X0.

python-control/control/timeresp.py

Lines 940 to 944 in 55730b6

     if isinstance(sys, TransferFunction) and np.any(X0 != 0):  
   warnings.warn(  
   "Non-zero initial condition given for transfer function system. "  
   "Internal conversion to state space used; may not be consistent "  
   "with given X0.")  

 

`force_response() takes the initial condition as initial state of the system. This is not necessarily equivalent to initial output of the transfer function.

p…

View full answer

bnavigator · 2023-01-05T21:09:28Z

bnavigator
Jan 5, 2023
Collaborator

The code you posted should have given you a UserWarning:

UserWarning: Non-zero initial condition given for transfer function system. Internal conversion to state space used; may not be consistent with given X0.

python-control/control/timeresp.py

Lines 940 to 944 in 55730b6

    
           if isinstance(sys, TransferFunction) and np.any(X0 != 0): 
        
               warnings.warn( 
        
                   "Non-zero initial condition given for transfer function system. " 
        
                   "Internal conversion to state space used; may not be consistent " 
        
                   "with given X0.")

`force_response() takes the initial condition as initial state of the system. This is not necessarily equivalent to initial output of the transfer function.

python-control/control/timeresp.py

Line 1007 in 55730b6

xout[:, 0] = X0

python-control/control/timeresp.py

Lines 1029 to 1051 in 55730b6

    
           # Algorithm: to integrate from time 0 to time dt, with linear 
        
               # interpolation between inputs u(0) = u0 and u(dt) = u1, we solve 
        
               #   xdot = A x + B u,        x(0) = x0 
        
               #   udot = (u1 - u0) / dt,   u(0) = u0. 
        
               # 
        
               # Solution is 
        
               #   [ x(dt) ]       [ A*dt  B*dt  0 ] [  x0   ] 
        
               #   [ u(dt) ] = exp [  0     0    I ] [  u0   ] 
        
               #   [u1 - u0]       [  0     0    0 ] [u1 - u0] 
        
               M = np.block([[A * dt, B * dt, np.zeros((n_states, n_inputs))], 
        
                            [np.zeros((n_inputs, n_states + n_inputs)), 
        
                             np.identity(n_inputs)], 
        
                            [np.zeros((n_inputs, n_states + 2 * n_inputs))]]) 
        
               expM = sp.linalg.expm(M) 
        
               Ad = expM[:n_states, :n_states] 
        
               Bd1 = expM[:n_states, n_states+n_inputs:] 
        
               Bd0 = expM[:n_states, n_states:n_states + n_inputs] - Bd1 
        
               for i in range(1, n_steps): 
        
                   xout[:, i] = (Ad @ xout[:, i-1] 
        
                                 + Bd0 @ U[:, i-1] + Bd1 @ U[:, i]) 
        
               yout = C @ xout + D @ U

2 replies

bnavigator Jan 5, 2023
Collaborator

gain=2
tc=20
bias=13.5

def fo_model(gain, tc):
         num=np.array([gain])
         den=np.array([tc,1])
         H = control.ss(control.tf(num,den))
         X0 = (pv[0]-bias)/H.C
         (t_ ,y) = control.forced_response(H, t, cv, X0=X0)
         return y+bias     
  
plt.plot(t, pv, 'b', label='PV')
plt.plot(t, cv, 'r', label='CV')
plt.plot(t, fo_model(gain,tc), 'g', label='FO')
plt.plot([0], pv[0], '+', label='initial condition')
plt.legend()
plt.show()

TheWicklowWolf Jan 6, 2023
Author

Exactly what I was looking for, excellent!

I have one follow up question (initial conditions when using a delay) below.
Thank you for your help, it is much appreciated.

TheWicklowWolf · 2023-01-06T12:29:39Z

TheWicklowWolf
Jan 6, 2023
Author

How do you correctly set the Initial conditions when using a delay (Padé approximation)?

t=[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
    18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
    36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
    54, 55, 56, 57, 58, 59]

cv=[10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20]

pv=[33.55,33.57,33.39,33.45,33.54,33.35,33.36,33.59,33.57,33.43,33.47,33.56
    ,33.35,34.4,35.4,36.26,37.13,37.95,38.54,39.23,39.94,40.77,41.31,41.92
    ,42.37,42.97,43.47,44.02,44.38,44.96,45.4,45.76,46.09,46.45,46.9,47.12
    ,47.36,47.65,47.98,48.18,48.5,48.79,48.94,49.19,49.43,49.78,49.94,50.04
    ,50.21,50.38,50.56,50.75,50.83,51.05,51.1,51.18,51.26,51.33,51.41,51.68]

import numpy as np
import matplotlib.pyplot as plt
import control

gain=2
tc=20
dt=5
bias=13.5

def fopdt_model(gain, tc):
    num=np.array([gain])
    den=np.array([tc,1])
    #Add Delay
    n_pade = 10
    (num_pade,den_pade) = control.pade(dt,n_pade)
    H_pade=control.tf(num_pade,den_pade)

    H_noDelay = control.ss(control.tf(num,den))
    X0 = (pv[0]-bias)/H_noDelay.C

    X0pade=np.append(np.zeros(n_pade),X0)
    H = control.series(H_pade,H_noDelay)
    (t_ ,y) = control.forced_response(H, t, cv, X0=X0pade)
    return y+bias  

plt.plot(t, pv, 'b', label='PV')
plt.plot(t, cv, 'r', label='CV')
plt.plot(t, fopdt_model(gain,tc), 'g', label='FOPDT')
plt.legend()
plt.show()

This is the result (I think it just requires the correct Initial Padé condition):

3 replies

bnavigator Jan 6, 2023
Collaborator

Totally depends on what you consider "correct" initial conditions.

TheWicklowWolf Jan 6, 2023
Author

I wasn't very clear, but I was hoping I could get the Padé approximation to be somewhat consistent with the Initial PV in the Plot above.
Maybe in on the wrong path but I think I need to set the initial Padé approximation to be none zero to get the desired result...

TheWicklowWolf Jan 6, 2023
Author

This is what I have so far:

def fopdt_model(gain, tc,dt):
    num=np.array([gain])
    den=np.array([tc,1])
    H_noDelay = control.ss(control.tf(num,den))
    X0 = (pv[0]-bias)/H_noDelay.C

    #Add Delay
    n_pade = 10
    (num_pade,den_pade) = control.pade(dt,n_pade)
    H_pade=control.ss(control.tf(num_pade,den_pade))
    X0pade=np.divide((pv[0]-bias),H_pade.C,where=H_pade.C!=0) #prevent infinity when dividing by zero

    X0combined=np.append(X0pade,X0)
    H = control.series(H_pade,H_noDelay)
    (t_ ,y) = control.forced_response(H, t, cv, X0=X0combined)
    return y+bias

But its actually a worse response:

TheWicklowWolf · 2023-01-06T20:18:41Z

TheWicklowWolf
Jan 6, 2023
Author

The Padé approximation works perfectly when the initial state is zero:
I would like to emulate that when the initial state is not zero, I'm just not sure how I can do that?

0 replies

roryyorke · 2023-01-08T15:38:14Z

roryyorke
Jan 8, 2023

The simplest way is to treat the linear model as changes from the initial steady-state, and so subtract off cv[0] and pv[0] when solving the linear model, and adding them back afterwards (see below).

You could solve for steady x0 given u0 (cv[0]), which is solving A x0 = -B u0 for x0, but then you need to check that y0 = C x0. matches your pv[0]. It looks like you were doing something of the sort already with bias.

t=[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
    18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
    36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
    54, 55, 56, 57, 58, 59]

cv=[10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20]

pv=[33.55,33.57,33.39,33.45,33.54,33.35,33.36,33.59,33.57,33.43,33.47,33.56
    ,33.35,34.4,35.4,36.26,37.13,37.95,38.54,39.23,39.94,40.77,41.31,41.92
    ,42.37,42.97,43.47,44.02,44.38,44.96,45.4,45.76,46.09,46.45,46.9,47.12
    ,47.36,47.65,47.98,48.18,48.5,48.79,48.94,49.19,49.43,49.78,49.94,50.04
    ,50.21,50.38,50.56,50.75,50.83,51.05,51.1,51.18,51.26,51.33,51.41,51.68]

import numpy as np
import matplotlib.pyplot as plt
import control

gain=2
tc=20
dt=5

def fopdt_model(gain, tc):
    num=np.array([gain])
    den=np.array([tc,1])
    #Add Delay
    n_pade = 10
    (num_pade,den_pade) = control.pade(dt,n_pade)
    H_pade=control.ss(control.tf(num_pade,den_pade))

    H_noDelay = control.ss(control.tf(num,den))

    H = control.series(H_pade,H_noDelay)

    (t_ ,y) = control.forced_response(H, t, np.array(cv) - cv[0])
    return y + pv[0]

plt.plot(t, pv, 'b', label='PV')
plt.plot(t, cv, 'r', label='CV')
plt.plot(t, fopdt_model(gain,tc), 'g', label='FOPDT')
plt.legend()
plt.show()

1 reply

TheWicklowWolf Jan 9, 2023
Author

That makes it much a cleaner implementation, thank you both for your help. It is much appreciated.

The end (keeping the bias as an external parameter) result:

t=[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
    18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
    36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
    54, 55, 56, 57, 58, 59]

cv=[10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
    20, 20, 20, 20, 20, 20]

pv=[33.55,33.57,33.39,33.45,33.54,33.35,33.36,33.59,33.57,33.43,33.47,33.56
    ,33.35,34.4,35.4,36.26,37.13,37.95,38.54,39.23,39.94,40.77,41.31,41.92
    ,42.37,42.97,43.47,44.02,44.38,44.96,45.4,45.76,46.09,46.45,46.9,47.12
    ,47.36,47.65,47.98,48.18,48.5,48.79,48.94,49.19,49.43,49.78,49.94,50.04
    ,50.21,50.38,50.56,50.75,50.83,51.05,51.1,51.18,51.26,51.33,51.41,51.68]

import numpy as np
import matplotlib.pyplot as plt
import control

gain=2
tc=20
dt=5
bias=13.5

def fopdt_model(gain, tc, dt, bias):
    num=np.array([gain])
    den=np.array([tc,1])
    #Add Delay
    n_pade = 10
    (num_pade,den_pade) = control.pade(dt,n_pade)
    H_pade=control.ss(control.tf(num_pade,den_pade))

    H_noDelay = control.ss(control.tf(num,den))

    H = control.series(H_pade,H_noDelay)

    (t_ ,y) = control.forced_response(H, t, np.array(cv) - cv[0])
    return y + bias + (cv[0]*gain)

plt.plot(t, pv, 'b', label='PV')
plt.plot(t, cv, 'r', label='CV')
plt.plot(t, fopdt_model(gain,tc,dt,bias), 'g', label='FOPDT')
plt.legend()
plt.show()

	if isinstance(sys, TransferFunction) and np.any(X0 != 0):
	warnings.warn(
	"Non-zero initial condition given for transfer function system. "
	"Internal conversion to state space used; may not be consistent "
	"with given X0.")

Issue with Non-Zero starting conditions when using control.forced_response - Not working as expected? #843

Uh oh!

Uh oh!

TheWicklowWolf Jan 5, 2023

Replies: 4 comments · 6 replies

Uh oh!

bnavigator Jan 5, 2023 Collaborator

Uh oh!

Uh oh!

bnavigator Jan 5, 2023 Collaborator

Uh oh!

TheWicklowWolf Jan 6, 2023 Author

Uh oh!

TheWicklowWolf Jan 6, 2023 Author

Uh oh!

bnavigator Jan 6, 2023 Collaborator

Uh oh!

TheWicklowWolf Jan 6, 2023 Author

Uh oh!

TheWicklowWolf Jan 6, 2023 Author

Uh oh!

TheWicklowWolf Jan 6, 2023 Author

Uh oh!

roryyorke Jan 8, 2023

Uh oh!

TheWicklowWolf Jan 9, 2023 Author

TheWicklowWolf
Jan 5, 2023

Replies: 4 comments 6 replies

bnavigator
Jan 5, 2023
Collaborator

bnavigator Jan 5, 2023
Collaborator

TheWicklowWolf Jan 6, 2023
Author

TheWicklowWolf
Jan 6, 2023
Author

bnavigator Jan 6, 2023
Collaborator

TheWicklowWolf Jan 6, 2023
Author

TheWicklowWolf Jan 6, 2023
Author

TheWicklowWolf
Jan 6, 2023
Author

roryyorke
Jan 8, 2023

TheWicklowWolf Jan 9, 2023
Author