Lab 2 – Braitenberg Lab – answers — v2
Lab 2 — sample solutions¶
Copyright By PowCoder代写 加微信 powcoder
Introduction to Braitenberg Vehicles¶
Ethology is the study of animal behaviour, and thus part of biology. Just as we can use computational methods to investigate population dynamics, genomics, phylogenetics, etc., we can use computational methods as part of ethology. This is an area of computational biology that overlaps at times with the cognitive sciences, including psychology, neuroscience, philosophy of mind and artificial intelligence. It’s worth noting the difference in goals between these different fields. Often in AI, the primary goal is to create impressive or useful technology, whereas the other cognitive sciences tend to be motivated more directly by trying to understand how intelligence is accomplished in Nature.
A famous (and very easy to read) book entitled Vehicles[1] by shows how some basic forms of intelligent and adaptive behaviour can be accomplished using very basic forms of feedback, where…(i) what is sensed affects what action is taken, and (ii) the action that is taken affects what is sensed in the next moment.
This ongoing feedback cycle of (a) affects (b) affects (a) affects (b) is a powerful concept, but it can be difficult/impossible to imagine what happens when there is feedback in a system—for whatever reason, it’s not something that our brains are good at doing.
This is a place where (computational) modelling can help. Instead of philosopher’s arm-chair experiments, it is possible to “instantiate your ideas,” i.e. to build mathematical models of the ideas you have and to use those models to investigate the implications of the assumptions and other facets of the model you have developed. The famous philospher has written about this role of modelling[2], as have other authors that have inspired me[3].
Along these lines, in this lab, you will develop a simulation of one of the ‘vehicles’ that Braitenberg developed in the course of his book.
[1] Braitenberg, V. (1986). Vehicles: Experiments in synthetic psychology. MIT Press.
[2] Dennett, D. (1995). Artificial life as philosophy. Artificial Life. An overview, 291-2.
[3] Di Paolo, E. A., Noble, J., & Bullock, S. (2000). Simulation Models as Opaque Thought Experiments. _Seventh International Conference on Artificial Life,_497–506.
Euler’s method for integration¶
In this lab, we will be using ordinary differential equations (ODE’s) to model the movement of robots. First we will introduce a method for numeric integration known as Euler’s method.
An ODE can be used to describe how a rate changes with over time. However, most ODE’s do not have general solutions, so we will need to use numeric methods to approximate the solution.
Euler’s method is defined as follows:
Given an ODE $f(t, y)$, and integration interval $(a, b)$
The approximate solution $\hat{y}$ using $N$ iterations can be computed by
Step 1 Define initial conditions $t_0 = a$, $w_0 = y_0$, number of iterations $N$ and step size $h = \frac{b – a}{N}$
Step 2 Starting at index $i = 1$
Step 3 Calculate time $t_i$, and approximation $w_i$ where
$t_i = t_0 + h \cdot i$
$w_i = w_{i-1} + h \cdot f(t_{i-1}, w_{i-1})$
Step 4 Check loop condition. If $i < N$, then increment $i$ and repeat Step 3.
Otherwise if $i = N$, return your solution $\hat{y} = (w_0, w_1, ..., w_N)$
See here for details on the intuition behind Euler's method.
Here, we will implement Euler's method to solve the equation below
$f(t, y) = \frac{dy}{dt} = \frac{1}{t^2} - \frac{y}{t} - y^2$,
with initial condition $y(1) = -1$, and time interval $1 \leq t \leq 2$.
Fill out the skeleton functions below. The integrate() function should print out Euler's approximation for each step of the iteration.
import numpy as np
import matplotlib.pyplot as plt
def dydt(t, y):
"""Differential equation goes here"""
result = 1 / (t**2) - y / t - y**2
return result
def integrate(y0, a, b, N):
"""Implement Euler's method here"""
h = (b - a) / N
t0 = a # initial time at t0
wi = y0 # y0 solution at t0
print("i = 0, ti = %s, wi = %s" % (ti, wi))
while i <= N:
wi = wi + h * dydt(ti, wi)
ti = t0 + h * i
print("i = %s, ti = %s, wi = %s" % (i, ti, wi))
Check that your implementation is working correctly using $N = 10$, and the conditions
$y_0 = -1, a = 1, b = 2$ and $t_0 = a$.
Your method should produce the following output.
Step Time Euler's approximation
$i$ $t_i$ $w_i$
0 1.0 -1.0
1 1.1 -0.9000000
2 1.2 -0.8165372
3 1.3 -0.7457213
4 1.4 -0.6847965
5 1.5 -0.6317567
6 1.6 -0.5851068
7 1.7 -0.5437101
8 1.8 -0.5066872
9 1.9 -0.4733469
10 2.0 -0.4431388
a = 1.0 # lower time interval
b = 2.0 # upper time interval
integrate(y0, a, b, N)
i = 0, ti = 1.0, wi = -1
i = 1, ti = 1.1, wi = -0.9
i = 2, ti = 1.2, wi = -0.8165371900826447
i = 3, ti = 1.3, wi = -0.745721278076786
i = 4, ti = 1.4, wi = -0.6847965276646216
i = 5, ti = 1.5, wi = -0.6317567102411785
i = 6, ti = 1.6, wi = -0.5851068058741311
i = 7, ti = 1.7000000000000002, wi = -0.5437101279350207
i = 8, ti = 1.8, wi = -0.5066871734303049
i = 9, ti = 1.9, wi = -0.4733468776585252
i = 10, ti = 2.0, wi = -0.4431387797315269
Once you are satisfied that your implementation is correct.
We can copy and adapt your integrate function for the next tasks.
Braitenberg Vehicles¶
Instead of considering in accurate detail the way that specific organisms move, Braitenberg imagines a set of 'vehicles' or robots. These are, for the most part, very simple. We shall consider a vehicle that has two powered wheels, and two light sensors. We might imagine the robot to look something like this:
But our model will simplify this robot significantly. In fact our model will have only three dynamic variables:
$x$ and $y$ describe the position of the robot in 2D space
$\alpha$ describes the direction that the robot is facing
The robot has two separate motors; one for each wheel. The left motor velocity is described by the term $l_m$ and the right motor by $r_m$.
These motors affect $x$, $y$, and $\alpha$ according to the following differential equations.
\begin{align}
\frac{dx}{dt} & = k_{m} \cos(\alpha) (l_m + r_m) \\
\frac{dy}{dt} & = k_{m} \sin(\alpha) (l_m + r_m) \\
\frac{d\alpha}{dt} & = k_{m} \frac{r_m-l_m}{2 k_r}
\end{align}where
$k_{m} = 10.0$ is a parameter describing the maximum speed of the robot
$k_r = 0.5$ is the radius of the robot
Fill in the differential equations in the methods below dxdt(), dydt(), dadt(). In the code, we will use 'a' to represent alpha.
For now we'll imagine the motors lm(), rm() are fixed to constant values.
"""Returns the state of the left motor. This function will be updated/replaced later."""
return 1.0
"""Returns the state of the right motor. This function will be updated/replaced later."""
return 1.0
The arguments (x,y, and a) are the current state of the system. Each function $\frac{dx}{dt}$, $\frac{dy}{dt}$, $\frac{d \alpha}{dt}$ should return the change in the relevant variable given the current state of the system, as defined by the three equations above.
k_m = 10.0
def dxdt(x, y, a, time):
the_rate_of_change_in_x = k_m * np.cos(a) * (lm() + rm())
return the_rate_of_change_in_x
def dydt(x, y, a, time):
the_rate_of_change_in_y = k_m * np.sin(a) * (lm() + rm())
return the_rate_of_change_in_y
def dadt(x, y, a, time):
the_rate_of_change_in_a = k_m * (rm() - lm()) / (2 * k_r)
return the_rate_of_change_in_a
Adapt your Euler integration function to numerically integrate the system of three differential equations above.
The function takes as arguments the initial values of the three variables $x_0, y_0, a_0$ and the simulation duration (end time), assuming start time $t_0 = 0$.
For this task, set the number of iterations at $N = 1500$.
The function should return 3 vectors, containing the approximated values at each step $\hat{x} = (x_0, x_1, ..., x_N)$, $\hat{y} = (y_0, y_1, ..., y_N)$ and $\hat{\alpha} = (\alpha_0, \alpha_1, ..., \alpha_N)$.
def integrate_system(x0, y0, a0, endtime, N):
t0 = 0 # initial time
h = (endtime - t0) / N
xi = x0 # x0 initial solution at t0
yi = y0 # y0 initial solution at t0
ai = a0 # a0 initial solution at t0
results = np.zeros((3, N + 1))
results[:, 0] = [x0, y0, a0]
while i <= N:
xi_next = xi + h * dxdt(xi, yi, ai, ti) # calculate state using previous states
yi_next = yi + h * dydt(xi, yi, ai, ti)
ai_next = ai + h * dadt(xi, yi, ai, ti)
ti = t0 + h * i
# print("i = %s, ti = %s, xi = %s, yi = %s, ai = %s" % (i, ti, xi, yi, ai))
results[0, i] = xi_next # store trajectory
results[1, i] = yi_next
results[2, i] = ai_next
xi = xi_next # update previous states
yi = yi_next
ai = ai_next
return results
y = integrate_system(1.0, 1.0, 2, 10.0, 1500)
Plot out the $x$, $y$, and $\alpha$ values of your robot. This is the trajectory of your robot.
You can select any value for the end time duration. Try using different numbers for the motor rates (edit the return value of the lm and rm functions). Do you get the behaviour you expect from your robot?
Use plt.plot(x) to plot the trajectory of variable $x$.
See the documentation of the plot function here: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html
plt.plot(y[0, :], y[1, :])
[
Additional Tasks¶
Ok, so you now have a simple simulation of a robot driving around, but its motor rates are unchanging, and so the paths that it moves through are not that interesting. Taking inspiration from Braitenberg, your next task is to make it so that the motor states are not constant, but are instead driven by the light-sensors that are on the robot.
We will start by assuming that there is a single source of light and that it is always located at the origin $(x,y) = (0,0)$
The following code returns the excitation of a light sensor that is located at the x,y position specified in the arguments that is facing in the direction specified by the angle argument. Note that the angle argument is specified in radians, not degrees.
The sensor works such that when it is directly facing the light it is maximally excited (a value of 1). This excitation falls off when the sensor does not face the light directly. The distance between the light and the sensor has no affect upon its excitation. There are lots of other ways that I could have done this.
You do not need to edit the code below, but you will use it for additional tasks 1 and 2.
def sensorAt(x,y,angle) :
# a vector from the sensor’s position to the origin (which is where the light is located)
dsquared = (x**2+y**2)
# normalize a vector from sensor to light
mag = np.sqrt(dsquared)
s2l_x = -x/mag
s2l_y = -y/mag
# unit vector indicating direction of sensor
h_x = np.cos(angle)
h_y = np.sin(angle)
attenuation_factor = (s2l_x*h_x + s2l_y*h_y)
attenuation_factor = max(0.0, attenuation_factor)
return attenuation_factor
The following code calculates the position and angle of the sensors to be able to use the sensorAt function to calculate the excitation of the two sensors. Again, you don’t need to edit this code, but you may be interested to take a look to see if you can understand how it works.
SENSOR_OFFSET = 2.0*np.pi/6.0
def ls(x,y,a) :
“””Given the position and orientation of the robot,
this function returns the excitation of the left sensor.”””
sensor_a = a + SENSOR_OFFSET
sensor_x = x + np.cos(sensor_a)*k_r
sensor_y = y + np.sin(sensor_a)*k_r
return sensorAt(sensor_x,sensor_y,sensor_a)
def rs(x,y,a) :
“””Given the position and orientation of the robot,
this function returns the excitation of the right sensor.”””
sensor_a = a – SENSOR_OFFSET
sensor_x = x + np.cos(sensor_a)*k_r
sensor_y = y + np.sin(sensor_a)*k_r
return sensorAt(sensor_x,sensor_y,sensor_a)
Additional Task 1¶
Using these two functions (ls and rs), rewrite the lm and rm functions to match the following equations.
\begin{align}
l_m &= -0.025 + r_s \\
r_m &= -0.020 + l_s
\end{align}Here we are creating a version of Braitenberg’s “aggression” vehicle.
You will have to also update the calls to these functions in your dxdt,dydt and dadt functions above, so that they include the x,y,and a variables that are now necessary for calculating the motor values (as they are necessary for calculating the sensor values).
def lm(x,y,a):
return -0.025 + rs(x, y, a)
def rm(x,y,a):
return -0.020 + ls(x, y, a)
k_m = 10.0
def dxdt(x, y, a, time):
the_rate_of_change_in_x = k_m * np.cos(a) * (lm(x, y, a) + rm(x, y, a))
return the_rate_of_change_in_x
def dydt(x, y, a, time):
the_rate_of_change_in_y = k_m * np.sin(a) * (lm(x, y, a) + rm(x, y, a))
return the_rate_of_change_in_y
def dadt(x, y, a, time):
the_rate_of_change_in_a = k_m * (rm(x, y, a) – lm(x, y, a)) / (2 * k_r)
return the_rate_of_change_in_a
Additional Task 2¶
Write a loop that generates 10 trajectories, and plots each one in a different color on the same figure. Confirm that the robots move toward the light (given enough time) from any reasonably close location.
Set the duration of the simulation to $50.0$, number of iterations to $N= 5000$ (or step size $h = 0.01$).
You may like to use the plt.plot functions below:
import matplotlib.pyplot as plt
import numpy as np
X = np.arange(1, 5, 0.01) # example X trajectory
Y = np.arange(1, 5, 0.01) # example Y trajectory
plt.plot(X, Y, ‘.’, alpha=0.1) # marker type ‘.’ for trajectory
plt.plot(X[-1], Y[-1], ‘o’) # marker type ‘o’ for end of trajectory
Your output should look something like this:
Experiment with plotting the trajectory data in different ways to understand how/why the robots move toward the light.
As time allows, experiment with different parameters and ways of hooking up the sensors to the motors to see what kinds of basic behaviours can be accomplished using these extremely simple form of feedback control.
What does the constant SENSOR_OFFSET do? How does it change the robot’s behaviour when you change it?
Why are the constants in the $l_m$ and $r_m$ equations different? What happens when you set them to the same value? Why?
endtime = 50.0
for trial in range(10) :
x0 = np.random.randn() * 10.0
y0 = np.random.randn() * 10.0
a0 = np.random.rand() * np.pi * 2.0
trajectory = integrate_system(x0, y0, a0, endtime, N)
plt.plot(trajectory[0, :], trajectory[1, :],’.’, alpha=0.1)
plt.plot(trajectory[0, -1], trajectory[1, -1], ‘o’)
plt.plot(0, 0, ‘x’)
plt.xlim(-r, r)
plt.ylim(-r, r)
plt.show()
Sample solutions
The SENSOR_OFFSET controls the ‘eye angle’ of the motor. A number of possible behaviours are possible, depending on the setup conditions for the simulation. When the $l_m$ and $r_m$ are set to the default assymetric values, a large sensor offset value will cause the trajectory radius to be further away from the light source.
The constants in the $l_m$ and $r_m$ equations are different so the left and right motors are slightly asymmetrical. This will mean that when the motor starts off away from the light source, the asymmetry will prevent it from moving purely in a straight line.
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com