My approach does what you want. You can see it by calling lower_bound.(acceleration) and upper_bound.(acceleration).
The key is to note the i in:
-max_acceleration + accel_offset[i] <= acceleration[i=1:2, 1:T] <= max_acceleration + accel_offset[I]
It results in
-max_acceleration + 0.0 <= acceleration[1, 1:T] <= max_acceleration + 0.0
for i=1 and
-max_acceleration - accel_g <= acceleration[2, 1:T] <= max_acceleration - accel_g
for i=2.