I have been working for awhile on how to best do this in a general way in python-control. The challenge is that the controller outputs a zero-order-hold “staircase” function that is incompatible with the way variable-step numerical integrators work. (They jump back and forth). I have not yet figured out a way to force any of scipy’s integrators to pick only the time steps you specify, but I am not sure I have exhausted all options yet.

In the meantime one way to simulate such a system is to repeatedly simulate your dynamical system (eg tanks2) over short intervals ofts usinginput_output_response With the constant inputs produced by your discrete-time controller at each ts.

This is cumbersome and slow though. What I have been working on is a an alternative that entails simple fixed-step Euler integration, at a short dt for accuracy, of the continuous-time plant. This entails creating a discrete-time model of the plant for a short dt, and a function that can create anonlineariosystem that implements the staircase function itself from a discrete-time controller (which has a longer ts sample time). Then, these can be connected usinginterconnect and simulated usinginput_output_response.

I’ll post here if I get something working.

FWIW, here's one approach, in which one callssolve_ivp once for each sample interval; within a given interval, the control signal is constant. One could extend this by concatenating thesol.t andsol.y results (remembering to eliminate common points) to generate a continuous time solution. One can also write code to combine the dense solution interpolants, a bit like Matlab'sodextend.

This is fairly raw; I haven't thought at all about how to apply this to the IOSystem classes.

importnumpyasnpfromscipy.integrateimportsolve_ivpimportmatplotlib.pyplotasplt# non-linear DEdefde(t,y,u):return-y-y**3+u# controller; discrete or continuous timedefprop_controller(r,y,k):returnk* (r-y)# desired final timetf_des=0.3# sample perioddt=1/10# number of samplesnsamples=max(5,int(np.ceil(tf_des/dt)))k=5# gainr=5# setpointy0=-1# initial state# discrete-time state (=output)yd=np.empty(nsamples)# discrete-time inputud=np.empty(nsamples)# initializeyd[0]=y0ud[0]=prop_controller(r,yd[0],k)fornsampleinrange(1,nsamples):# construct new DE function for each sample intervalf=lambdat,y:de(t,y,ud[nsample-1])sol=solve_ivp(f, [nsample*dt, (nsample+1)*dt], [yd[nsample-1]])yd[nsample]=sol.y[0,-1]ud[nsample]=prop_controller(r,yd[nsample],k)# continuous-time controlled systemdefcontrol_de(t,y,r,k):returnde(t,y,prop_controller(r,y,k))# continuous-time solutioncont_sol=solve_ivp(lambdat,y:control_de(t,y,r,k),                     [0,nsamples*dt], [y0])# compare continous-time and discrete-timeplt.clf()plt.stairs(yd,np.arange(nsamples+1)*dt,baseline=None)plt.plot(cont_sol.t,cont_sol.y[0])

dt=10 result:

dt=100 result:

At the end, I modified the solution by@roryyorke: created a helper function that correctly integrates within sampling time and a wrapper for nlsys that evaluates the model in discrete-time. This will give correct results at sampling times and can be used without further modifications in existing code. If a more detailed solution is needed between sampling times, it is easy to integrate the system using the resulting piece-wise constant closed-loop input and solve_ivp with more data points.

def nonlinear_plant_dynamics(t, x, u, params):    F1, F2, k11, k22 = 0.5, 0.6, 0.8, 0.5    h1, h2 = x    q01, q02 = u    h2 = 0 if h2 < 0 else h2   # numerical issues    h1 = h2 if h1 < h2 else h1    dh1 = (q01 - k11 * math.sqrt(h1 - h2)) / F1    dh2 = (q02 + k11 * math.sqrt(h1 - h2) - k22 * math.sqrt(h2)) / F2    return np.array((dh1, dh2))def nonlinear_plant_output(t, x, u, params):    return x# simulate plant dynamics within sampling perioddef discrete_nonlinear_plant_dynamics(t, x, u, params):    f = lambda t, x: nonlinear_plant_dynamics(t, x, u, params)    sol = solve_ivp(f, [t, t+ts], x)    xn = sol.y[:, -1]    return nonlinear_plant_output(t, xn, u, params)# nonlinear_plant = ct.nlsys(nonlinear_plant_dynamics, nonlinear_plant_output,  inputs=['q01','q02'], outputs=['h1', 'h2'], states=['h1', 'h2'], name='tanks2i')discrete_nonlinear_plant = ct.nlsys(discrete_nonlinear_plant_dynamics, nonlinear_plant_output, dt=ts, inputs=['q01','q02'], outputs=['h1', 'h2'], states=['h1', 'h2'], name='tanks2i_d')

Yes, I studied the example in the docs and played originally with it as well. Discrete-time delays can indeed be modeled using state-space formulation. As you mention, the integration precision was less accurate and I wanted the exact results as I had in Simulink. In addition, this proposed solution requires only one small additional function and the whole code can remain the same - it integrates very well with the rest of the python-control package.