2
$\begingroup$

I want to numerically simulate the motion of a point particle (unit mass) in $N$-dimensional space according to the force $f(\vec{x}, t)$ and bounded by a set of linear constraints (equalities as well as inequalities):

$$A_i \vec{x} = \vec{c_i}$$

$$B_j \vec{x} \leq \vec{d_j}$$

The behavior of the particle when it comes in contact with a "wall" formed by an inequality is a perfectly inelastic collision with no friction. So when it hits a wall, it slides along said wall instead of bouncing or sticking.

What I'm thinking now is something like this:

for step in range(total_time//dt):
    p += f(x, step*dt) * dt
    x += p * dt
    if violate_constraints(x):
        x, p = correct(x, p)

where x and p are numpy arrays. It is easy to check which constraints $\vec{x}$ violates, but I'm not sure how to implement the correct(x, p) part.

Considering the "corner" cases where the particle comes in contact with 2 or more walls in a single step, I think the correct behavior for correct(x, p) is to put the position $\vec{x}$ back to the closest point in the "allowed" space, and project $\vec{p}$ onto the intersection of subspaces according to which constraints $\vec{x}$ (before correction) broke. However,

  1. I have no proof if this is correct;
  2. I don't know implement the correct(x, p) function, especially considering I want a working solution for arbitrary constraints.
$\endgroup$

1 Answer 1

4
$\begingroup$

The correction you are looking for is a projection $$ \newcommand{\norm}[1]{\left\lVert#1\right\rVert} \DeclareMathOperator*{\argmin}{arg\,min} P_{C} : \mathbb{R}^{d} \rightarrow C, \, \, \pmb{x} \mapsto \argmin_{\pmb{y} \in C} \norm{\pmb{x} - \pmb{y}}_{2} $$

onto a convex polytope $C$ given by $$ C = \lbrace \pmb{x} \in \mathbb{R}^{d} : A \pmb{x} = \pmb{c}, B \pmb{x} \leq \pmb{d} \rbrace $$

Then $P_{C}$ is determined by minimizing the following quadratic function subject to given linear constraints: \begin{align*} & \text{minimize} \ \frac{1}{2} \norm{\pmb{x} - \pmb{y}}_{2}^{2} \\ & \text{subject to: } \tag{P}\\ & \qquad \pmb{a}_{i}^{T} \pmb{x} = c_{i} \quad \forall i = 1, \dots, m \\ & \qquad \pmb{b}_{i}^{T} \pmb{x} \leq d_{j} \quad \forall j = 1, \dots, n \end{align*}

Optimization problem (P) belongs to the more general class of problems called Quadratic Programming (QP). It is thus very simple to implement $P_{C}$ by using off-the-shelf solver.

Here is a simple implementation using CVXPY in Python:

import cvxpy as cp
import numpy as np

def project(y, A, c):
    d = A.shape[1]
    x = cp.Variable(d)
    obj = cp.Minimize(cp.sum_squares(x - y))
    constr = [A @ x <= c]
    prob = cp.Problem(obj, constr)
    prob.solve()
    return np.array(x.value).squeeze()

Example. Given a convex polytope $C$ and a point $\pmb{x} = (1.77, 0.23)^T \not\in C$, as shown in Figure (a) below, then running:

A = np.array([[2, -1], [-2, 1], [0, 1], [-1, -1], [1, -1]])
c = np.array([3, 1, 2, -1, 1])
x = np.array([1.77, 0.23])
project(x, A, c)

we get the projection $P_{C}(\pmb{x}) \in C$ of $\pmb{x}$ onto $C$:

array([1.5, 0.5])

Finally, to simulate the motion of a particle as you specified in your question. Say that particle starts at position $\pmb{x}_{0} = (1, 0.5)^T$ and a constant force $f$ acts on it so that it moves by some vector $\pmb{v} = (0.4, 0.1)^T$ at each time step:

def f(x):
    v = np.array([0.4, 0.1])
    return x + v

Then, by running:

nsteps = 8
x0 = np.array([1, 0.5])
xs = np.zeros((nsteps + 1, 2))
xs[0, ] = x0
for i in range(nsteps):
    y = f(xs[i])
    xs[i + 1, ] = project(y, A, c)

we collect in the list xs the positions of the particle shown in Figure (b) below:

particle projection

Remark. You can try to adapt the code for this simple example for your specific problem. It should work for dimensions higher than 2. Also, note that each equality constraint can be replaced by a pair of inequalities: $$ ax + by = c \iff ax + by \geq c \text{ and } ax + by \leq c $$ so the function project can also be used for this case. By adding each equality constraint, we are intersecting with affine linear space and get a polytope $C$ of dimension $d-1$, one less than the dimension of the ambient space $\mathbb{R}^d$.

Update. Regarding your first point of considering "corner" cases:

  • When $d > 2$ we should also consider the edge case, literally. If particle hits the edge and leaves the set $C$, then by convexity of $C$, the particle is projected back by $P_{C}$ to the same edge. This is true, if we use small enough step size (two inequality constraints will be active). Otherwise, we can overshoot and particle will be projected onto the corner of $C$ (at least 3 inequality constraints will be active).
  • As for the corner case, this is easier to see. The particle will be projected back to the corner, by convexity of $C$.
  • All in all, you should observe the described behaviour when simulating the motion of a particle in $C$ using projection operator $P_{C}$.
  • You should only use small enough step size, otherwise the particle will not slide along the walls, but will be projected straight into the corners of $C$.

To find on which boundaries the particle is: We just need to check which constraints are active after applying the projection operator.

This is easy, since the boundaries, or we can call them walls, are linear equations determined by matrix $A$ and vector $c$.

Using the above example, say that the particle is moved by force $f$ to position $\pmb{x} = (1.77, 0.23)^T$ and projected back onto $C$ by $P_C$ to position $\pmb{y} = (1.5, 0.5)^T$.

Then we can check which constraints are active by running:

A @ y == c

to get:

array([False, False, False, False,  True])

Or for a nicer result:

[bnd[0] + 1 for bnd in np.transpose(np.nonzero(A @ y == c))]

to get:

[5]

Meaning that the particle is on the 5th boundary $x - y = 1$, as we can also observe in Figure (a) above.

$\endgroup$
2
  • 1
    $\begingroup$ Thanks for the thorough answer! I'm not familiar with cvxpy, is there any easy way to know which boundaries x is on after it hits a wall/walls? $\endgroup$
    – haoyu
    Sep 16, 2022 at 12:06
  • 1
    $\begingroup$ Yes, I added an update how to find which boundaries the particle is on. $\endgroup$ Sep 16, 2022 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.