|
| 1 | +""" |
| 2 | +Mix Integer Optimization based path planner |
| 3 | +
|
| 4 | +author: Atsushi Sakai |
| 5 | +""" |
| 6 | + |
| 7 | + |
| 8 | +import cvxpy |
| 9 | +import math |
| 10 | +import numpy as np |
| 11 | +import matplotlib.pyplot as plt |
| 12 | +from matplotrecorder import matplotrecorder |
| 13 | +matplotrecorder.donothing = True |
| 14 | + |
| 15 | + |
| 16 | +# parameter |
| 17 | +A = np.matrix([[1.0, 0.0], |
| 18 | + [0.0, 1.0]]) |
| 19 | +B = np.matrix([[1.0, 1.0], |
| 20 | + [0.0, 1.0]]) |
| 21 | +q = np.matrix([[1.0], |
| 22 | + [1.0]]) |
| 23 | +r = np.matrix([[0.1], |
| 24 | + [0.1]]) |
| 25 | + |
| 26 | +u_max = 0.1 # input constraint |
| 27 | +T = 30 # horizon length |
| 28 | +M = 10000.0 # weight parameter of obstacle avoidanse |
| 29 | + |
| 30 | + |
| 31 | +def control(s1, gs, ob): |
| 32 | + |
| 33 | + w = cvxpy.Variable(2, T) |
| 34 | + v = cvxpy.Variable(2, T) |
| 35 | + s = cvxpy.Variable(2, T) |
| 36 | + u = cvxpy.Variable(2, T) |
| 37 | + nob = len(ob) |
| 38 | + o = cvxpy.Bool(4 * nob, T) |
| 39 | + |
| 40 | + constraints = [cvxpy.abs(u) <= u_max] |
| 41 | + constraints.append(s[:, 0] == s1) |
| 42 | + |
| 43 | + obj = [] |
| 44 | + for t in range(T): |
| 45 | + constraints.append(s[:, t] - gs <= w[:, t]) |
| 46 | + constraints.append(-s[:, t] + gs <= w[:, t]) |
| 47 | + constraints.append(u[:, t] <= v[:, t]) |
| 48 | + constraints.append(-u[:, t] <= v[:, t]) |
| 49 | + |
| 50 | + obj.append(t * q.T * w[:, t] + r.T * v[:, t]) |
| 51 | + |
| 52 | + # obstable avoidanse |
| 53 | + for io in range(nob): |
| 54 | + ind = io * 4 |
| 55 | + constraints.append(sum(o[ind:ind + 4, t]) <= 3) |
| 56 | + constraints.append(s[0, t] <= ob[io, 0] + M * o[ind + 0, t]) |
| 57 | + constraints.append(-s[0, t] <= -ob[io, 1] + M * o[ind + 1, t]) |
| 58 | + constraints.append(s[1, t] <= ob[io, 2] + M * o[ind + 2, t]) |
| 59 | + constraints.append(-s[1, t] <= -ob[io, 3] + M * o[ind + 3, t]) |
| 60 | + |
| 61 | + for t in range(T - 1): |
| 62 | + constraints.append(s[:, t + 1] == A * s[:, t] + B * u[:, t]) |
| 63 | + |
| 64 | + objective = cvxpy.Minimize(sum(obj)) |
| 65 | + |
| 66 | + prob = cvxpy.Problem(objective, constraints) |
| 67 | + |
| 68 | + prob.solve(solver=cvxpy.GUROBI) |
| 69 | + |
| 70 | + s_p = s.value |
| 71 | + u_p = u.value |
| 72 | + print("status:" + prob.status) |
| 73 | + |
| 74 | + return s_p, u_p |
| 75 | + |
| 76 | + |
| 77 | +def plot_obstacle(ob): |
| 78 | + for i in range(len(ob)): |
| 79 | + x = [ob[i, 0], ob[i, 1], ob[i, 1], ob[i, 0], ob[i, 0]] |
| 80 | + y = [ob[i, 2], ob[i, 2], ob[i, 3], ob[i, 3], ob[i, 2]] |
| 81 | + plt.plot(x, y, "-g") |
| 82 | + |
| 83 | + |
| 84 | +def main(): |
| 85 | + print(__file__ + " start!!") |
| 86 | + |
| 87 | + s = np.matrix([10.0, 5.0]).T # init state |
| 88 | + gs = np.matrix([5.0, 7.0]).T # goal state |
| 89 | + |
| 90 | + ob = np.matrix([[7.0, 8.0, 3.0, 8.0], |
| 91 | + [5.5, 6.0, 6.0, 10.0]]) # [xmin xmax ymin ymax] |
| 92 | + # ob = np.matrix([[7.0, 8.0, 3.0, 8.0]]) |
| 93 | + |
| 94 | + h_sx = [] |
| 95 | + h_sy = [] |
| 96 | + |
| 97 | + for i in range(10000): |
| 98 | + print("time:", i) |
| 99 | + s_p, u_p = control(s, gs, ob) |
| 100 | + |
| 101 | + s = A * s + B * u_p[:, 0] # simulation |
| 102 | + |
| 103 | + if(math.sqrt((gs[0] - s[0]) ** 2 + (gs[1] - s[1]) ** 2) <= 0.1): |
| 104 | + print("Goal!!!") |
| 105 | + break |
| 106 | + |
| 107 | + h_sx.append(s[0, 0]) |
| 108 | + h_sy.append(s[1, 0]) |
| 109 | + |
| 110 | + plt.cla() |
| 111 | + plt.plot(gs[0], gs[1], "*r") |
| 112 | + plot_obstacle(ob) |
| 113 | + plt.plot(s_p[0, :], s_p[1, :], "xb") |
| 114 | + plt.plot(h_sx, h_sy, "-b") |
| 115 | + plt.plot(s[0], s[1], "or") |
| 116 | + plt.axis("equal") |
| 117 | + plt.grid(True) |
| 118 | + plt.pause(0.0001) |
| 119 | + matplotrecorder.save_frame() |
| 120 | + |
| 121 | + matplotrecorder.save_movie("animation.gif", 0.1) |
| 122 | + |
| 123 | + |
| 124 | +if __name__ == '__main__': |
| 125 | + main() |
0 commit comments