Code snippets used in ACC 2025 tool paper¶
This python notebook goes over the various code snippets and examples demonstrated in the ACC 2025 toolpaper describing pycvxset.
Vinod, A.P., "pycvxset: A Python package for convex set manipulation", American Control Conference (ACC), 2025.
License information: pycvxset code is released under AGPL-3.0-or-later license, as found in the LICENSE.md file. The documentation for pycvxset is released under CC-BY-4.0 license, as found in the LICENSES/CC-BY-4.0.md.
This notebook follows the same organization as the ACC 2025 toolpaper describing pycvxset:
The links below are for the HTML page in the documentation website, and they will not work in the Jupyter notebook.
- Code block 1: Introduction to polytopes
- Code block 2: Introduction to constrained zonotopes
- Code block 3: Introduction to ellipsoids
- Figure 1 and Code block 4: Visualizing polytopes and polytopic approximations
- Figure 2: Spreading points on a 3d unit sphere
- Figure 3 and Code block 5: Projection of the point
- Figure 4: Centering for a polytope and a constrained zonotope
- Figure 5 and Code block 6: orthogonal projection of a unit l1 norm ball
- Code block 7: Checking for equality
- Code block 8: Reachability analysis using pycvxset
- Figure 6: Robust controllable set computation for a double integrator
- Figure 7: Robust controllable set computation for a spacecraft rendezvous example
Code block 1: Introduction to polytopes¶
The following code snippet creates a polytope in V-Rep and $3$-dimensional simplex in H-Rep, prints the description of the polytope along with its vertices.
import numpy as np
from pycvxset import Polytope
V = [[-1, 0.5], [-1, 1], [1, 1], [1, -1], [0.5, -1]]
P1 = Polytope(V=V)
print("P1 is a", repr(P1))
A, b = -np.eye(3), np.zeros((3,))
Ae, be = [1, 1, 1], 1
P2 = Polytope(A=A, b=b, Ae=Ae, be=be)
print("P2 is a", P2)
print("Vertices of P2 are:\n", P2.V)
print("P2 is a", P2)
P1 is a Polytope in R^2 in only V-Rep In V-rep: 5 vertices P2 is a Polytope in R^3 in only H-Rep Vertices of P2 are: [[0. 0. 1.] [0. 1. 0.] [1. 0. 0.]] P2 is a Polytope in R^3 in H-Rep and V-Rep
The call $\texttt{P2.V}$ in Line 11 triggers a vertex enumeration internally as seen from the $\texttt{print}$ statements for $\texttt{P2}$.
Code block 2: Introduction to constrained zonotopes¶
The following code snippet creates a constrained zonotope from the polytope defined before as well as a box.
from pycvxset import ConstrainedZonotope
C1 = ConstrainedZonotope(polytope=P1)
print("C1 is a", repr(C1))
print("P1 is a", repr(P1))
C2 = ConstrainedZonotope(lb=[-1, -1], ub=[1, 1])
print("C2 is a", repr(C2))
C1 is a Constrained Zonotope in R^2 with latent dimension 5 and 1 equality constraint P1 is a Polytope in R^2 in only V-Rep In V-rep: 5 vertices C2 is a Constrained Zonotope in R^2 that is a zonotope with latent dimension 2
Note that $\texttt{pycvxset}$ detects that $\texttt{C2}$ is a zonotope.
Code block 3: Introduction to ellipsoids¶
$\texttt{pycvxset}$ supports ellipsoids that are full-dimensional or degenerate. The following code snippet creates two ellipsoids of different forms.
from pycvxset import Ellipsoid
E1 = Ellipsoid(c=[2, -1], Q=np.diag([1, 4]))
print("E1 is an", E1)
E2 = Ellipsoid(c=[0, 1, 0], G=np.diag([1, 2, 3]))
print("E2 is an", E2)
E1 is an Ellipsoid in R^2 E2 is an Ellipsoid in R^3
Figure 1 and Code block 4: Visualizing polytopes and polytopic approximations¶
We plot 2D and 3D polytopes and plot polytopic approximations of constrained zonotopes and ellipsoids.
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = [5, 3]
plt.rcParams["figure.dpi"] = 150
%matplotlib widget
plt.figure()
ax = plt.subplot(131, projection="3d")
P2.plot(ax=ax) # Plot polytope
ax.view_init(elev=30, azim=-15)
ax.set_aspect("equal")
ax.set_title("Polytope")
ax = plt.subplot(132) # Plot const. zonotope
C1.plot(ax=ax, vertex_args={"visible": True})
ax.set_aspect("equal")
ax.set_title("Constrained Zonotope")
ax = plt.subplot(133) # Plot ellipsoid
E1.plot(ax=ax, patch_args={"facecolor": "pink"})
ax.set_aspect("equal")
ax.set_title("Ellipsoid")
plt.subplots_adjust(wspace=0.5)
Figure 2: Spreading points on a 3D unit sphere¶
We spread points on a 3D unit sphere. We use these points for generating (inner- and outer-) polytopic approximations within pycvxset. These approximations are useful for plotting sets that do not have a direct conversion to V-Rep.
from pycvxset import spread_points_on_a_unit_sphere
n = 3
D = 20
n_points = 2 * n + (2**n) * D
dir_vectors = spread_points_on_a_unit_sphere(n_dim=n, n_points=n_points)[0]
plt.figure()
ax = plt.subplot(111, projection="3d")
ax.scatter(dir_vectors[:, 0], dir_vectors[:, 1], dir_vectors[:, 2])
ax.set_aspect("equal")
ax.set_xticks([-1, 0, 1])
ax.set_yticks([-1, 0, 1])
ax.set_zticks([-1, 0, 1])
ax.set_xticklabels([-1, 0, 1], fontsize=20)
ax.set_yticklabels([-1, 0, 1], fontsize=20)
ax.set_zticklabels([-1, 0, 1], fontsize=20)
ax.set_xlabel("x", fontsize=20, labelpad=20)
ax.set_ylabel("y", fontsize=20, labelpad=20)
ax.set_zlabel("z", fontsize=20, labelpad=20);
Figure 3 and Code block 5: Projection of the point¶
We compute the Euclidean projection and the distance of a point $[1, 1, 1]$ on the polytope $\texttt{P2}$.
projection, d = P2.project(x=[1, 1, 1], p=2)
# Plotting
point = np.array([1, 1, 1])
patch_args = {"label": "Polytope"}
ax = P2.plot(patch_args=patch_args)[0]
ax.view_init(elev=15, azim=-45)
ax.scatter(*point, color="red", label="Point $x$")
ax.scatter(*projection[0], label="Projection")
ax.plot(
[point[0], projection[0, 0]],
[point[1], projection[0, 1]],
[point[2], projection[0, 2]],
"k--",
label=f"Dist.: {d[0]:1.2f}",
)
ax.legend(bbox_to_anchor=(1.35, 0.8), fontsize=15)
ax.set_xticks([0, 0.5, 1])
ax.set_yticks([0, 0.5, 1])
ax.set_zticks([0, 0.5, 1])
ax.set_xticklabels([0, 0.5, 1], fontsize=13)
ax.set_yticklabels([0, 0.5, 1], fontsize=13)
ax.set_zticklabels([0, 0.5, 1], fontsize=13)
ax.set_xlabel("x", fontsize=13)
ax.set_ylabel("y", fontsize=13)
ax.set_zlabel("z", fontsize=13)
ax.set_aspect("equal")
plt.subplots_adjust(right=0.4);
Figure 4: Centering for a polytope and a constrained zonotope¶
We illustrate centering and bounding sets for the polytope $\texttt{P1}$ and the constrained zonotope $\texttt{C1}$.
from pycvxset import is_polytope
for C1_type, plot_sets in zip(["vrep", "hrep"], [False, True]):
print(f"Using CZ obtained from {C1_type:s} polytope")
if plot_sets:
plt.figure()
if C1_type == "hrep":
C1 = ConstrainedZonotope(polytope=Polytope(A=P1.A, b=P1.b))
elif C1_type == "vrep":
C1 = ConstrainedZonotope(polytope=Polytope(V=P1.V))
for index, cset in enumerate([P1, C1]):
rect = Polytope.deflate_rectangle(cset)
cheby_ball = Ellipsoid.inflate_ball(cset)
max_vol_ell = Ellipsoid.inflate(cset)
patch_args = {
"facecolor": "pink",
"linewidth": 3,
"linestyle": ":",
"label": "Bounding rectangle",
}
if plot_sets:
ax = plt.subplot(121 + index)
rect.plot(ax, patch_args=patch_args)[0]
cset.plot(ax=ax, patch_args={"label": "Set"})[0]
cheby_ball.plot(
ax=ax,
center_args={"color": "k", "label": "Chebyshev. center"},
patch_args={
"facecolor": "gray",
"label": "Chebyshev ball",
},
)
max_vol_ell.plot(
ax=ax,
center_args={"color": "gold", "label": "Ellipsoid center"},
patch_args={
"facecolor": "gold",
"alpha": 0.6,
"label": "Max. vol. ellipsoid",
},
)
if index == 1:
handles, labels = ax.get_legend_handles_labels()
handles = [handles[1], handles[0], *handles[2:]]
labels = [labels[1], labels[0], *labels[2:]]
ax.legend(handles, labels, bbox_to_anchor=(1, 1.02))
ax.set_aspect("equal")
if is_polytope(cset):
ax.set_title("Polytope")
else:
ax.set_title("Constrained Zonotope")
print(
f"Set {type(cset).__name__}: Chebyshev radius (in each dim.): "
f"{np.array2string(np.diag(cheby_ball.G), formatter={'float_kind':'{:1.2f}'.format}):s}"
)
print(
f"Set {type(cset).__name__}: MVIE volume: {max_vol_ell.volume():1.2f}"
)
plt.subplots_adjust(right=0.6, wspace=0.3)
Using CZ obtained from vrep polytope Set Polytope: Chebyshev radius (in each dim.): [0.73 0.73] Set Polytope: MVIE volume: 1.89 Set ConstrainedZonotope: Chebyshev radius (in each dim.): [0.61 0.61] Set ConstrainedZonotope: MVIE volume: 1.22 Using CZ obtained from hrep polytope
/home/vinod/workspace/git/PaperCodes/general/pycvxset/pycvxset/ConstrainedZonotope/operations_unary.py:78: UserWarning: This function returns a sub-optimal but feasible solution for Chebyshev centering.
warnings.warn("This function returns a sub-optimal but feasible solution for Chebyshev centering.", UserWarning)
/home/vinod/workspace/git/PaperCodes/general/pycvxset/pycvxset/ConstrainedZonotope/operations_unary.py:146: UserWarning: This function returns a sub-optimal but feasible solution for maximum volume inscribed ellipsoid-based centering.
warnings.warn(
Set Polytope: Chebyshev radius (in each dim.): [0.73 0.73] Set Polytope: MVIE volume: 1.89
Set ConstrainedZonotope: Chebyshev radius (in each dim.): [0.73 0.73] Set ConstrainedZonotope: MVIE volume: 1.89
Figure 5 and Code block 6: Orthogonal projection of a unit l1-norm ball¶
In the following code snippet, we compute the projection of a $3$-dimensional unit $\ell_1$-norm ball.
V = np.vstack((np.eye(3), -np.eye(3)))
l1ball = Polytope(V=V)
ball2D = l1ball.projection(project_away_dims=2)
# Plotting the sets
fig = plt.figure()
ax = fig.add_subplot(1, 2, 1, projection="3d")
ax.view_init(**{"elev": 18, "azim": -100})
l1ball.plot(ax=ax, patch_args={"facecolor": "lightblue"})
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_aspect("equal")
ax.set_title("3D polytope")
ax2d = fig.add_subplot(1, 2, 2)
ball2D.plot(ax=ax2d)
ax2d.set_xlabel("x")
ax2d.set_ylabel("y")
ax2d.set_aspect("equal")
ax2d.set_title("Projection on x-y plane")
ax2d.grid()
plt.subplots_adjust(wspace=0.4)
Code block 7: Checking for equality¶
The following code block checks for equality between polytopes and constrained zonotopes. All equality checks are performed by executing a pair of containment checks. Note that polytope containing another polytope is a convex program. Similarly, constrained zonotope containing a V-Rep polytope is also a convex program. However, constrained zonotope containing another constrained zonotope or a H-Rep polytope is a bilinear program (and requires GUROBI).
print("C1 is a", C1)
print("P1 is a", P1)
print("Are C1 and P1 equal?", C1 == P1)
lb, ub, p, q = [-1, -1], [1, 1], [-1, -1], 0.5
P1a = Polytope(lb=lb, ub=ub).intersection_with_halfspaces(A=p, b=q)
C1a = ConstrainedZonotope(lb=lb, ub=ub).intersection_with_halfspaces(
A=p, b=q
)
print("C1a is a", C1a)
print("P1a is a", P1a)
print("Are P1 and P1a equal?", P1a == P1a)
try:
print("Are C1 and C1a equal?", C1 == C1a)
except ValueError as err:
# Error handling if GUROBI is not set up properly
print("Skipped due to error ", print(err))
C1 is a Constrained Zonotope in R^2 P1 is a Polytope in R^2 in H-Rep and V-Rep Are C1 and P1 equal? True C1a is a Constrained Zonotope in R^2 P1a is a Polytope in R^2 in only H-Rep Are P1 and P1a equal? True
Set parameter ServerPassword
Set parameter TokenServer to value "atlas3.merl.com"
Are C1 and C1a equal? True
Code block 8: Reachability analysis using pycvxset¶
Consider a discrete-time linear time-invariant system with additive uncertainty,
\begin{align}
x_{t+1} = A x_t + B u_t + F w_t,
\end{align}
with state $x_t\in\mathbb{R}^n$, input $u_t\in\mathcal{U}\subset\mathbb{R}^m$, disturbance $w_t\in\mathcal{W}\subset\mathbb{R}^p$, and appropriate matrices $A,B,F$.
We assume that the input set $\mathcal{U}$ and disturbance set $\mathcal{W}$ are convex and compact sets.
Given a horizon $N\in\mathbb{N}$, a polytopic safe set $\mathcal{S}\subset\mathbb{R}^n$ and a polytopic target set $\mathcal{T}\subset\mathbb{R}^n$, a \emph{$N$-step robust controllable set} is the set of initial states that can be robustly driven, through a time-varying control law, to the target set in $N$ steps, while satisfying input and state constraints for all possible disturbances.
Formally, we define the $N$-step RC set as $\mathcal{K}_0$ via the following set recursion for $t\in\{0, 1, \ldots, N-1\}$:
\begin{align}
\mathcal{K}_t = \mathcal{S} \cap \left({A^{-1}{((\mathcal{K}_{t+1} \ominus F\mathcal{W}) \oplus (-B \mathcal{U}))}}\right),
\end{align}
with $\mathcal{K}_N\triangleq \mathcal{T}$. We implement the recursion in pycvxset with the following Python function \texttt{get_rcs}.
def get_rcs(S_U, S_W, S_S, S_T, A, B, F, N):
S_K = [None] * (N + 1)
S_K[-1], S_FW, S_BU = S_T, F @ S_W, (-B) @ S_U
for t in range(N - 1, -1, -1):
S_temp = (S_K[t + 1] - S_FW) + S_BU
S_K[t] = S_S.intersection(S_temp @ A)
return S_K
Figure 6: Robust controllable set computation for a double integrator¶
We use get_rcs function defined above to compute the robust controllable set for a double integrator. Observe that the
RC set computed using constrained zonotope is slightly smaller than the set computed using polytopes, due to the
inner-approximation used in Pontryagin difference
import time
# Dynamics
N = 30
sampling_time = 0.1
A = np.array([[1, sampling_time], [0, 1]])
B = np.array([[(sampling_time**2) / 2], [sampling_time]])
F = B
disturbance_scaling = 0.4
# Define the terminal set and input set (Polytope)
plot_skip = N // 3
set_definitions_list = [
[
Polytope(lb=[-0.25, -0.1], ub=[0.25, 0.1]), # S_T
Polytope(lb=[-1.0, -0.5], ub=[1.0, 0.5]), # S_S
Polytope(lb=-1, ub=1), # S_U
],
[
ConstrainedZonotope(lb=[-0.25, -0.1], ub=[0.25, 0.1]), # S_T
ConstrainedZonotope(lb=[-1.0, -0.5], ub=[1.0, 0.5]), # S_S
ConstrainedZonotope(lb=-1, ub=1), # S_U
],
]
# Compute the sets
elapsed_time_compute = [None] * 2
rcs = [None] * 2
for index, (S_T, S_S, S_U) in enumerate(set_definitions_list):
S_W = disturbance_scaling * S_U
start_time = time.time()
rcs[index] = get_rcs(S_U, S_W, S_S, S_T, A, B, F, N)
elapsed_time_compute[index] = time.time() - start_time
# Plot the sets
elapsed_time_plot = [None] * 2
plt.figure()
for index, (S_K, (S_T, S_S, S_U)) in enumerate(
zip(rcs, set_definitions_list)
):
start_time = time.time()
ax = plt.subplot(121 + index)
S_S.plot(
ax=ax,
patch_args={
"alpha": 1,
"facecolor": "yellow",
"label": "Safe set $\mathcal{S}$",
},
)
start_time = time.time()
S_K[0].plot(
ax=ax,
patch_args={
"alpha": 1,
"facecolor": "tab:blue",
"label": "$\mathcal{K}_0$",
},
)
for t, fc in zip(
range(plot_skip, N, plot_skip),
["lightseagreen", "palegreen", "aquamarine"],
):
S_K[t].plot(
ax=ax,
patch_args={
"alpha": 1,
"facecolor": fc,
"label": "$\mathcal{K}_{" + f"{t:d}" + "}$",
},
)
elapsed_time_plot[index] = time.time() - start_time
S_T.plot(
ax=ax,
patch_args={
"alpha": 1,
"facecolor": "white",
"label": f"Target set $\mathcal{{K}}_{{{N}}}=\mathcal{{T}}$",
},
)
ax.set_aspect("equal")
ax.set_xlabel("Position")
if index == 0:
ax.set_ylabel("Velocity")
else:
handles, labels = ax.get_legend_handles_labels()
new_ordering = [0, 4, 3, 2, 1]
new_handle = [handles[i] for i in new_ordering]
new_label = [labels[i] for i in new_ordering]
ax.legend(new_handle, new_label, bbox_to_anchor=(1.25, 1.5))
if is_polytope(S_S):
ax.set_title("Polytope")
else:
ax.set_title("Const. Zono.")
elapsed_time_plot[index] = time.time() - start_time
plt.subplots_adjust(right=0.6, wspace=0.4)
plt.tight_layout()
print(
f"Time taken for to compute Polytope: {elapsed_time_compute[0]:1.3f} s"
)
print(f"Time taken for to plot Polytope: {elapsed_time_plot[0]:1.3f} s")
print(
f"Time taken for to compute ConstrainedZonotope: {elapsed_time_compute[1]:1.3f} s"
)
print(
f"Time taken for to plot ConstrainedZonotope: {elapsed_time_plot[1]:1.3f} s"
)
plt.show()
/home/vinod/workspace/git/PaperCodes/general/pycvxset/pycvxset/ConstrainedZonotope/operations_binary.py:644: UserWarning: This function returns an inner-approximation of the Pontryagin difference with a zonotope. warnings.warn(
Time taken for to compute Polytope: 1.887 s Time taken for to plot Polytope: 0.149 s Time taken for to compute ConstrainedZonotope: 0.397 s Time taken for to plot ConstrainedZonotope: 3.583 s
Figure 7: Robust controllable set computation for a spacecraft rendezvous example¶
We now demonstrate a practical application of pycvxset where ConstrainedZonotope class provides scalability and
numerical stability over Polytope class for the computation of RC set. It also uses an ellipsoidal uncertainty set
defined using Ellipsoid class.
Specifically, we consider the problem of spacecraft rendezvous.
# Dynamics
system_A = np.array(
[
[1.0014, 0, 29.9953, 0.9225],
[0, 1.0000, -0.9225, 29.9811],
[0, 0, 0.9995, 0.0615],
[0, 0, -0.0615, 0.9981],
]
)
system_B = np.array(
[
[1.4999, 0.0308],
[-0.0308, 1.4995],
[0.1000, 0.0031],
[-0.0031, 0.0999],
]
)
system_F = np.eye(4)
los_max_y = 1 * (10**3) # Units m
los_max_velocity = 5e-1 # Units m/s
line_of_sight_cone_A = [
[1, 1, 0, 0],
[-1, 1, 0, 0],
[0, -1, 0, 0],
[0, 0, 1, 0],
[0, 0, -1, 0],
[0, 0, 0, 1],
[0, 0, 0, -1],
]
line_of_sight_cone_b = [
0,
0,
los_max_y,
los_max_velocity,
los_max_velocity,
los_max_velocity,
los_max_velocity,
]
line_of_sight_cone = Polytope(
A=line_of_sight_cone_A, b=line_of_sight_cone_b
)
line_of_sight_cone_constrained_zonotope = ConstrainedZonotope(
polytope=line_of_sight_cone
)
target_position_error = 0.2 * (10**3) # Units m
target_velocity_error = 0.1 # Units m/s
target_set = Polytope(
lb=[
-target_position_error,
-target_position_error,
-target_velocity_error,
-target_velocity_error,
],
ub=[
target_position_error,
0,
target_velocity_error,
target_velocity_error,
],
)
target_set_constrained_zonotope = ConstrainedZonotope(
lb=[
-target_position_error,
-target_position_error,
-target_velocity_error,
-target_velocity_error,
],
ub=[
target_position_error,
0,
target_velocity_error,
target_velocity_error,
],
)
max_acceleration = 0.2
input_constraints = Polytope(
lb=[-max_acceleration, -max_acceleration],
ub=[max_acceleration, max_acceleration],
)
input_constraints_constrained_zonotope = ConstrainedZonotope(
lb=[-max_acceleration, -max_acceleration],
ub=[max_acceleration, max_acceleration],
)
# disturbance_set = ConstrainedZonotope(c=[0,0,0,0], G=None)
disturbance_set = Ellipsoid(
c=[0, 0, 0, 0], G=np.diag([1e-2, 1e-2, 1e-4, 1e-4])
)
N = 50
We now compute the robust controllable sets.
start_time = time.time()
brs_list = get_rcs(
input_constraints_constrained_zonotope,
disturbance_set,
line_of_sight_cone_constrained_zonotope,
target_set_constrained_zonotope,
system_A,
system_B,
system_F,
N,
)
print(
f"Elapsed time for constrained zonotope computation: {time.time() - start_time:1.2f} seconds"
)
/home/vinod/workspace/git/PaperCodes/general/pycvxset/pycvxset/ConstrainedZonotope/operations_binary.py:636: UserWarning: This function returns an inner-approximation of the Pontryagin difference with an ellipsoid. warnings.warn(
Elapsed time for constrained zonotope computation: 24.59 seconds
We now plot the sets. For ease of plotting (and accuracy), we consider 3D projection of the polytope counterparts of the various sets, and then plot the robust controllable set (in constrained zonotope form) within them.
target_set_3D = [None] * 2
line_of_sight_cone_3D = [None] * 2
brs_3D = [None] * 2
for index, v_index in enumerate([2, 3]):
target_set_3D[index] = target_set.projection(v_index)
line_of_sight_cone_3D[index] = line_of_sight_cone.projection(v_index)
brs_3D[index] = (
brs_list[0]
.slice(v_index, 0)
.projection(v_index)
.polytopic_inner_approximation()
)
Generating the plot.
plt.figure()
for index, v_index in enumerate([2, 3]):
ax = plt.subplot(121 + index, projection="3d")
line_of_sight_cone_3D[index].plot(
ax=ax,
patch_args={
"facecolor": "yellow",
"label": "Safe set",
"alpha": 0.2,
},
)
brs_3D[index].plot(
ax=ax,
patch_args={
"facecolor": "tab:blue",
"label": "RC set",
"alpha": 1,
"edgecolor": "k",
},
)
target_set_3D[index].plot(
ax=ax,
patch_args={
"facecolor": "lightgreen",
"label": "Target set",
"alpha": 0.5,
},
)
ax.grid()
if v_index == 2:
# ax.set_title('$v_x=0$')
ax.set_zlabel("$v_y$ (m/s)")
[handles, labels] = ax.get_legend_handles_labels()
new_ordering = [0, 2, 1]
new_handles = [handles[i] for i in new_ordering]
new_labels = [labels[i] for i in new_ordering]
ax.legend(
new_handles,
new_labels,
loc="lower center",
ncols=3,
bbox_to_anchor=(1.25, -0.45),
)
ax.view_init(elev=20, azim=130)
else:
# ax.set_title('$v_y=0$')
ax.set_zlabel("$v_x$ (m/s)")
ax.view_init(elev=20, azim=130)
xticks_list = [-los_max_y, 0, los_max_y]
ax.set_xticks(xticks_list, [f"{int(v/1e3):d}" for v in xticks_list])
ax.set_xlabel("x (km)")
yticks_list = [-los_max_y, 0]
ax.set_yticks(yticks_list, [f"{int(v/1e3):d}" for v in yticks_list])
ax.set_ylabel("y (km)")
plt.subplots_adjust(wspace=0.45, right=0.85, bottom=0.2, left=0.05)