Squared Distance Problem

In this example, we demonstrate how to formulate and solve a general optimization problem within the OPTIMake framework, i.e., a problem with the number of stages set to 1.

When the number of stages is set to 1, the OPTIMake modeling interface reduces to:

$$\begin{split} &\quad \quad \quad \min_{v_1} l(v_1, p) \\ &\begin{split} \text{subject to} &\quad f_s(v_1, p) = 0,\\ &\quad g(v_1, p) \geq 0 \end{split} \end{split}$$

info

• The above optimization problem is a general optimization problem formulation, where the objective function, equality constraints, and inequality constraints are defined through the modeling interface functions objective, start_equality, and inequality, respectively.
• When the number of stages is 1, the multi-stage modeling interface functions cannot be called: end_objective, end_equality, equality.

warning

OPTIMake is not designed for solving large-scale general optimization problems. OPTIMake can efficiently model and solve general optimization problems (single-stage problems) when the following conditions are met:

• Variable/constraint dimensions are less than 100
• Fixed dimensions and sparsity structure

Problem Description

The minimum distance between a circle and a rectangle (when the circle and rectangle intersect, the minimum distance is 0) can be formulated as an optimization problem.

Let the center of the rectangle be $(x,y)$, with length and width $l,w$ respectively, the circle centered at the origin with radius $r$.

Suppose $(x_r,y_r)$ and $(x_o,y_o)$ are points inside the rectangle and the circle, respectively. Then the minimum distance problem can be described as minimizing the distance between $(x_r,y_r)$ and $(x_o,y_o)$, i.e.,

$$\begin{split} &\quad \quad \quad d_{\min}^2 =\min (x_r - x_o)^2 + (y_r - y_o)^2 \\ &\begin{split} \text{subject to} &\quad (x_r,y_r) \in \text{Rectangle},\\ &\quad (x_o,y_o) \in \text{Circle}\\ \end{split} \end{split}$$

The constraint that $(x_o,y_o)$ lies inside the circle can be described by the following quadratic constraint:

$$x_o^2 + y_o^2 \leq r^2$$

The constraint that $(x_r,y_r)$ lies inside the rectangle can be described in several ways. Here we use the vertices of the rectangle. Let the four vertices of the rectangle be $(v_1,v_2,v_3,v_4)$. Then a point inside the rectangle can be described as a weighted average of the 4 vertices, i.e.,

$$\begin{split} (x_r,y_r) = &\theta_1 v_1 + \theta_2 v_2 + \theta_3 v_3 + \theta_4 v_4 \\ &\theta_1 + \theta_2 +\theta_3 + \theta_4= 1\\ &\theta_1, \theta_2, \theta_3, \theta_4 \geq 0 \end{split}$$

Modeling

The above optimization problem is modeled in OPTIMake as follows:

prob = multi_stage_problem('squared_distance', 1)

# rectangle parameters
# center: (x, y), rotation: phi
# length: l, width: w
l, w = prob.parameters(['l', 'w'], stage_dependent=False)
x, y, phi = prob.parameters(['x', 'y', 'phi'], stage_dependent=False)
# circle parameters
# center: (0, 0), radius: r
r = prob.parameters(['r'], stage_dependent=False)

# point inside the circle
xo = prob.variable('xo')
yo = prob.variable('yo')
theta1 = prob.variable('theta1', hard_lowerbound=0.0)
theta2 = prob.variable('theta2', hard_lowerbound=0.0)
theta3 = prob.variable('theta3', hard_lowerbound=0.0)
theta4 = prob.variable('theta4', hard_lowerbound=0.0)

# rectangle vertices
V = Matrix([[+l / 2.0, +w / 2.0],
            [+l / 2.0, -w / 2.0],
            [-l / 2.0, -w / 2.0],
            [-l / 2.0, +w / 2.0]])
# rotation matrix
R = Matrix([[cos(phi), -sin(phi)],
            [sin(phi),  cos(phi)]])
# rotated rectangle vertices
V = V * R
theta = Matrix([theta1, theta2, theta3, theta4])

# point inside the rectangle
pos_rec = V.transpose() * theta + Matrix([[x], [y]])
xr, yr = pos_rec[0], pos_rec[1]

obj = general_objective((xr - xo)**2 + (yr - yo)**2)
prob.objective(obj)

seq = general_equality([theta1 + theta2 + theta3 + theta4 - 1.0])
prob.start_equality(seq)

ineq = general_inequality(
    expr = [xo**2 + yo**2 - r**2],
    sign = ['<='],
    bound = [0.0])
prob.inequality(ineq)

option = codegen_option()
option.default_tolerance_level = 'high'
codegen = code_generator()
codegen.codegen(prob, option)

Solving

#include "squared_distance_prob.h"
#include "squared_distance_solver.h"
#include <stdio.h>

int main(void)
{
    Squared_distance_Problem prob;
    Squared_distance_Option option;
    Squared_distance_WorkSpace ws;
    Squared_distance_Output output;
    squared_distance_init(&prob, &option, &ws);

    /* params */
    prob.param[SQUARED_DISTANCE_PARAM_X] = 2.0; /* x */
    prob.param[SQUARED_DISTANCE_PARAM_Y] = 2.0; /* y */
    prob.param[SQUARED_DISTANCE_PARAM_PHI] = 0.2; /* phi */
    prob.param[SQUARED_DISTANCE_PARAM_L] = 4.0; /* l */
    prob.param[SQUARED_DISTANCE_PARAM_W] = 2.0; /* w */
    prob.param[SQUARED_DISTANCE_PARAM_R] = 1.0; /* r */

    /* option */
    option.print_level = 1;
    int solve_status = squared_distance_solve(&prob, &option, &ws, &output);

    printf("solve_status = %d\n", solve_status);
    printf("squared distance = %f\n", output.obj);
    printf("(xo, yo) = (%f, %f)\n", \
        output.primal.var[0][SQUARED_DISTANCE_VAR_XO], output.primal.var[0][SQUARED_DISTANCE_VAR_YO]);

    return 0;
}

Results

When $x = 2$, $y = 2$, $\phi = 0.2$, $l = 4$, $w = 2$, $r = 1$, the computed $d_{\min}^2 = 0.127786$ (can be obtained via output.obj):

When $x = 2$, $y = 2$, $\phi = 0.0$, $l = 2$, $w = 2$, $r = 1$, the computed $d_{\min}^2 = 0.171573$:

When $x = 1.5$, $y = 1.5$, $\phi = 0.0$, $l = 2$, $w = 2$, $r = 1$, i.e., when the circle and rectangle intersect, the computed $d_{\min}^2 = 0$: