""" Ported by SCL, 27 Oct 2012. ...from combination of predprey.m and predprey_rh.m The original top-of-file comments appear in function docstrings. """ import numpy as np def predprey_rh(x, t, rnom, predprey_K, predprey_xd, predprey_ud): """ predprey_rh.m - feedback function for predator prey system RMM, 13 Sep 06 This function computes a simple linear control law, with saturation, for the predator prey system """ # Compute the feedback control law rh = (-np.dot(predprey_K, (x - predprey_xd)) + predprey_ud)[0] # Saturate the control input between 0 and 4*rnom if rh < 0: rh = 0 if rh > 4*rnom: rh = 4*rnom return rh def predprey(x, t, r = 1.6, d = 0.56, b = 0.6, k = 125, a = 3.2, c = 50, control_args=None): """ predprey.m - updated predator prey model for CDS 101/110 RMM, Oct 07 This file contains a simple model for a predator prey system. The model is developed following the treatment in the book by Murray (the other one). The states for these equations are x(1) # of rabbits (prey) x(2) # of foxes (predator) The equations allow for a quadratic term in both the predator and prey species that cause the population to fall off if you get too many of a given species. These terms are turned off by default. """ # Extract the states x1 = x[0] x2 = x[1] # Check to see if a feedback law has been defined (when r < 0) if r < 0 and (control_args is not None): # Assume that the user has defined a function that we can call # Add this value to the nominal value of the parameter (renegated) r = -r + predprey_rh(x, t, -r, *control_args) dx1 = r * x1 * (1 - x1/k) - a * x2 * x1/(c + x1) dx2 = b *a *x2 * x1 / (c + x1) - d * x2 return np.array([dx1, dx2])