Cloud Controls & Simulation Toolbox¶
The Cloud Controls & Simulation Toolbox (CCST) is a browser-based environment for controls engineering, simulation, and orbital mechanics. The backend is built with Python (FastAPI), using NumPy and SciPy for mathematical computations. The frontend terminal uses jQuery Terminal with Bokeh for interactive plotting. Simulations can run locally or as distributed ROS2 nodes in Docker or Kubernetes, with optional 3D visualization via Foxglove Studio.
All commands use arrow syntax to name outputs: command args -> output_name. Use gdisplay to view any stored variable (matrices, systems, plots, 3D scenes). Scripts (.ccst files) can be run with the run command.
Contact
For questions and bugs please email: greg.hayrapetyan AT gmail.com
Matrices and Display¶
Define a matrix using MATLAB-style or JSON syntax:
Display variables with gdisplay A (graphical LaTeX rendering) or display A (console text).
Calculate eigenvectors and eigenvalues:
Set display precision (0-16 decimal places, or full):
Running Scripts¶
Execute a .ccst script file (lines starting with # are comments):
Controls System Commands¶
Let us analyze longitudinal dynamics representative of a transport aircraft, trimmed at \(V_0 = 250\) ft/s and flying at a low altitude. The example is taken from 'Robust and Adaptive Control: With Aerospace Applications' by Eugene Lavretsky and Kevin Wise.
Defining the System¶
The linearized dynamics are given by the matrix:
To enter this matrix in CCST:
matrix '[[-0.038, 18.984, 0, -32.174], [-0.001,-0.632, 1, 0], [0, -0.759, -0.518, 0], [0, 0, 1, 0]]' -> A
View the matrix with gdisplay A.
Calculate eigenvectors and eigenvalues:
Control Matrices¶
Define the input, output, and feedthrough matrices:
matrix '[[10.1, 0], [0, -0.0086], [0.025, -0.011], [0,0]]' -> B
matrix '[[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]' -> C
matrix '[[0, 0], [0, 0], [0, 0], [0, 0]]' -> D
Creating State-Space System¶
Define the state-space system (\(\dot{x} = Ax + Bu\), \(y = Cx + Du\)):
View the system with gdisplay G. Generate a block diagram with:
Labeling States, Inputs, and Outputs¶
Label the second input as "elevator" and the second output as "alpha":
States can also be labeled:
Simulation¶
Simulate the step response to a 0.5 radian elevator deflection for 5 seconds:
Simulate response to initial conditions with a constant input:
Closed-Loop Control¶
Define a controller and arrange into a feedback loop:
PID Control¶
Wrap a PID controller around a plant output:
where the arguments are: plant, output name, \(K_p\), \(K_i\), \(K_d\).Transfer Functions and Poles¶
Compute the transfer function and plot poles:
Add zeros to also show zeros: poles G zeros -> pz_plot.
LQR, Kalman, and LQG¶
LQR State Feedback¶
Design an optimal state feedback gain by solving the algebraic Riccati equation. Define weighting matrices \(Q\) (state cost) and \(R\) (control cost):
Form the closed-loop system with state feedback (\(u = r - Kx\)):
Check closed-loop stability:
Kalman Filter¶
Define a stochastic system with process and measurement noise:
whereG_noise is the noise input matrix, Qn is process noise covariance, Rn is measurement noise covariance.
Design the optimal Kalman filter gain:
Create an observer (state estimator):
LQG Controller¶
Combine LQR state feedback and Kalman observer into an LQG controller:
System Augmentation¶
Add an integrator on an output (for zero steady-state error):
Add a tracking error output (\(e = r - y\)):
Add saturation (actuator limits):
Expose an internal signal as a plottable output:
Create a linear combination of outputs:
Cascade Control¶
Build inner/outer loop cascade architecture:
General Simulation¶
The sim command simulates any system (linear or nonlinear) with flexible input options.
Basic Simulation¶
Simulate a system from initial conditions for a given duration:
Specify a constant input and timestep:
Plot specific outputs by name:
Nonlinear Dynamics¶
Define arbitrary nonlinear systems using Python/NumPy expressions:
This defines a damped pendulum: \(\dot{\theta} = \omega\), \(\dot{\omega} = -\sin(\theta) - 0.1\omega\).
Simulate it:
Systems with inputs:
dynamics "[x[1], -np.sin(x[0]) - 0.1*x[1] + u[0]]" 2 1 -> forced_pendulum
sim forced_pendulum x0:[0,0] t:10 u:[0.5] -> forced_sim
Scheduled Inputs¶
Use target and chain to create piecewise input schedules (see Rendezvous & Docking tab), then pass them to sim:
3D Visualization¶
Create an interactive 3D scene from simulation results:
The scene includes a trajectory trail, play/pause controls, and orbit camera.
Rendezvous and Docking¶
CCST includes a complete set of tools for spacecraft proximity operations using Clohessy-Wiltshire (CW) relative motion dynamics.
Setup: CW Plant¶
Define the linearized relative motion dynamics in the LVLH frame for a 400 km LEO orbit (\(n = 0.00113\) rad/s):
This creates cw_plant (6-state linear system), an LQR controller K_rdv, and the closed-loop system rdv_cl.
Closed-Loop Rendezvous¶
Simulate a V-bar approach from 1 km behind and 100 m below:
sim rdv_cl x0:[-0.1,-1.0,0,0,0,0] t:3000 u:[0,0,0] dt:1 x y z -> vbar_pos
gdisplay vbar_pos
scene3d vbar_pos -> vbar_3d
gdisplay vbar_3d
Two-Impulse Targeting¶
Solve for the impulsive burns needed to transfer between two states:
For finite-duration burns (e.g., 30-second burns):
For thrust-limited burns:
Multi-Phase Docking¶
Chain multiple targeting phases with hold points:
target cw_plant x0:[-0.1,-1.0,0,0,0,0] xf:[0,-0.2,0,0,0,0] t:1400 duration:30 -> phase1
target cw_plant x0:[0,-0.2,0,0,0,0] xf:[0,-0.02,0,0,0,0] t:800 duration:30 -> phase2
target cw_plant x0:[0,-0.02,0,0,0,0] xf:[0,-0.006,0,0,0,0] t:400 duration:20 -> phase3
chain phase1 hold:300 phase2 hold:200 phase3 -> docking_sched
Simulate the complete approach:
sim cw_plant x0:[-0.1,-1.0,0,0,0,0] t:3100 schedule:docking_sched dt:1 x y z -> docking_traj
gdisplay docking_traj
scene3d docking_traj -> docking_3d
gdisplay docking_3d
Run the full docking scenario:
Available RDV Scripts¶
| Script | Scenario |
|---|---|
rdv_setup.ccst |
CW plant + LQR controller setup |
rdv_vbar.ccst |
V-bar approach with LQR |
rdv_outofplane.ccst |
Out-of-plane correction |
rdv_ellipse.ccst |
Safety ellipse departure |
rdv_targeting.ccst |
Two-impulse open-loop targeting |
rdv_finitethrust.ccst |
Finite-thrust maneuver |
rdv_attitude.ccst |
Attitude control |
rdv_coupled.ccst |
Coupled translation + attitude cascade |
rdv_docking.ccst |
Multi-phase docking with finite burns |
rendezvous.ccst |
Run all scenarios sequentially |
ROS2 Distributed Simulation¶
The ros2_sim command runs a simulation as distributed ROS2 nodes, with each block in the system's block diagram running as a separate node. This supports real-time pacing and 3D visualization via Foxglove Studio.
Basic ROS2 Simulation¶
Foxglove Visualization¶
Add the foxglove flag to publish visualization markers:
ros2_sim cw_plant x0:[-0.1,-1.0,0,0,0,0] t:3100 schedule:docking_sched dt:1 foxglove rtf:10 scale:100 x y z -> docking_ros2
Parameters:
foxglove: enable Foxglove bridge onws://localhost:8765rtf:<factor>: real-time factor (0 = fast as possible, 1 = real-time, 10 = 10x speed)scale:<factor>: position scaling for visualization (default 1000, i.e. km to m)frame:<var>:<parent>:<child>: add TF frame transforms
Connect Foxglove Studio to ws://localhost:8765 to view the 3D visualization in real time.
Managing Simulations¶
List running simulations:
Stop all running simulations:
Kubernetes Deployment¶
When deployed to Kubernetes (CCST_SIM_BACKEND=k8s), simulations run as K8s Jobs with S3-based data transfer. The system automatically handles job creation, data upload/download, and cleanup. See the K8s Deployment Guide for setup instructions.
Path Planning Toolbox¶
For speed the toolbox uses Python code for some of the functionality and C++ code for others.
Local Planner and 2D Car Dynamics¶
Algorithms like Dijkstra's or RRT/RDT are often very useful as global path planners, however there are situations when dynamic constraints have to be taken into account. There are several ways to model these dynamics. One could describe the car with the following state variables: the position of its center of mass in the Cartesian coordinates, \(x\), \(y\), its angle of orientation, \(\xi\), its velocity, \(v\), and its mass, \(m\). The inputs that can be controlled will be the amount of thrust applied, \(\delta T\), the amount of braking, \(\delta B\), and the angle by which the front wheels are turned, \(\delta w\).
Deriving the Dynamics¶
The velocity components in the \(x\) and \(y\) directions can be simply obtained as vector projections of \(v\) onto the coordinate axis using trigonometry:
We'll also assume that the mass changes as a function of velocity and the throttle, i.e.:
This can be specified more explicitly depending on the properties of the actual car.
Next we use Newton's Second law, \(F = m a = m \dot{v}\), written in terms of the velocity \(v\). To keep things a bit general we will assume that the total thrust provided by the car is \(T(v, \delta T, \delta B)\), a function of the current velocity, as well as throttle and brake amount. We also assume that the component of \(T\) perpendicular to the tire has an equal friction force in the opposite direction to cancel it, leaving the parallel component, \(T \cos(\delta w)\) as the net force. In turn, this force has components \(T \cos(\delta w) \cos(\delta w)\) in the direction of forward acceleration, \(\dot{v}\) and \(T \cos(\delta w) \sin(\delta w)\) in the tangential direction. It follows that:
where the last term is the drag due to air resistance, often approximated with \(C_D \rho v^2 S/2\), where \(C_D\) is the drag coefficient, \(\rho\) is the density of air, \(S\) is the area of the surface. We'll keep this as \(D(v)\) to keep things general.
The effect of the tangential force \(T \cos(\delta w) \sin(\delta w)\) is slightly more technical to derive using Newton's Law. To keep things relatively simple, we'll use a slightly different argument to obtain an expression for the angular velocity \(\dot{\xi}\). Assume that in time \(\Delta t\) the base of the car moved \(v \Delta t\) in original direction, \(\xi\), while the orientation \(\xi\) changed by \(\Delta \xi\).
Then, using the law of sines from trigonometry yields:
or after simplification:
Note that by the mean value theorem \(\sin(\delta w) - \sin(\delta w - \Delta \xi) = \cos(\delta w_*) \Delta \xi\) for some value of \(\delta w_*\) between \(\delta w - \Delta \xi\) and \(\delta w\). Substituting the expression for \(\sin(\delta w - \Delta \xi)\) results in:
or after simplification:
Finally, taking the limit as \(\Delta t \rightarrow 0\) (and consequently \(\delta w_* \rightarrow \delta w\)), we have:
Complete System of Equations¶
Let's summarize our derived system of differential equations:
Minimum Jerk Trajectories¶
To calculate a 5 second minimum jerk lane-change trajectory that takes a car traveling at 5m/s to the lane 4 meters to the left at the position 20 meters ahead we can run:
Here the arguments specify initial position \([x,y]\), initial velocity \([v_x, v_y]\), initial acceleration \([a_x, a_y]\), final position, final velocity, final acceleration, and the time to perform the maneuver.
Multiple Trajectory Generation¶
In real world trajectory generation oftentimes it is not enough to simply generate one optimal trajectory. In practice, one way to approach motion planning for lane changes or turns is to generate many possible trajectories by varying the final conditions and then selecting the best trajectory based on a more complicated cost function that could include, for example, distance to collisions, distance to the center of the lane, fuel consumption, etc.
Running the command:
generates 20 trajectories by varying the final velocities and positions, and chooses the best path based on ensuring forward direction of motion and jerk minimization. The additional parameters above specify number of grid points (1000), a flag that specifies that many trajectories should be sampled (True), and the number of trajectories to sample.
Implementation
My implementation of these solvers can be found here: C++ Code
The BVP solver is generic enough to handle not just the Euler-Lagrange equations corresponding to jerk minimization, but analogous functionals with arbitrary high derivatives and arbitrary number of waypoints.
For a more general collocation based BVP solver see: C++ Code
Optimal Stopping Control¶
To calculate the optimal control history to bring a car traveling at 20 m/s to a complete stop in 500 m we can run:
Dijkstra's Algorithm¶
One way to use Dijkstra's algorithm is to represent the map as a grid using an \(m\) by \(n\) matrix. To avoid entering a large matrix by hand the toolbox allows to upload any image whose grayscale can be thresholded to produce the obstacles.
Workflow¶
-
Upload Image: Expand the menu at the top right corner, above the console screen. Give a name to this image, for example
map1. This will be the handle for working with the image from the console. -
Display Image: After uploading the file we can display the uploaded image with
gdisplay map1. -
Generate Grid: To generate a grid by thresholding grayscale values of the image we can use:
wheremap1is the handle of the uploaded image and the last two parameters specify the height and width of the grid. The matrix representing the grid will be saved undermap1_gridname with values at each cell representing the grayscale value of the image. -
Create Obstacles: Entering
This sets cells with value larger than 2 to 1 (representing obstacle region) and cells with values less than or equal to 2 to 0 (representing free space). We can display the resulting matrix withgdisplay map1_griddisplays this matrix. To turn these values into free space and obstacle regions needed to run Dijkstra's algorithm we can enter:gdisplay map2_grid. -
Run Dijkstra's: Finally, running the following finds a path from the (11,0) position on the grid to the (8,6) position:
Orbital Mechanics Toolbox¶
The toolbox uses C++ code under the hood for performance.
Lambert's Problem¶
Consider a satellite to be located at \(r_1 = 5000 \hat{I} + 10000 \hat{J} + 2100 \hat{K}\) km. We want to design a trajectory that will put the satellite at \(r_2 = -14600 \hat{I} + 2500 \hat{J} + 7000 \hat{K}\) km after one hour. This is achieved by solving Lambert's problem, by calling:
Here we specified the two positions, time, whether we want the prograde trajectory, \(\mu\), and a name to save the results under. To see the orbital elements associated with the calculated orbit we simply type:
Optimal Circular Orbit Transfer¶
Consider establishing the optimal transfer trajectory to get from a given circular orbit around the earth to the largest radial position in a given amount of time using constant thrust and the thrust angle as the control. This problem can be posed as a classic optimal control problem of finding extrema of a functional of the form:
where the running cost \(L \equiv 0\) and the terminal cost \(K := r(t_f)\) is the final radial position. Precisely:
subject to the classic two-body equations of motion in polar coordinates:
where \(r\) is the radial position of the vehicle, \(u\) and \(v\) are radial and tangential velocities, \(T\) is the thrust applied at the angle of \(\alpha\) radians. In addition, \(m_0\) is the initial vehicle mass and \(\dot{m}\) is the mass flow rate.
Reference
See James M. Longuski, José J. Guzmán, John E. Prussing, Optimal Control with Aerospace Applications, Springer-Verlag New York (2014)
Boundary Value Problem Formulation¶
By the use of Pontryagin's Maximum Principle, this problem can be reduced to solving a boundary value problem:
for the rescaled variables:
the co-states \(p_r, p_{\theta}, p_u, p_v\) and subject to the boundary conditions:
Numerical Solution¶
To solve the problem with:
- \(\mu = 3.986004418 \times 10^5\) km³/s²
- Initial mass \(m_0 = 1500\) kg
- Specific impulse \(I_{sp} = 6000\) s
- Thrust \(T = 20\) N
- Initial orbit radius \(6678\) km
- Final time \(t_f = 0.5\) days
- Using 5 hour timesteps and 500 grid points for each timestep
we run:
The above solves the problem using a three stage Lobatto collocation scheme:
Similarly, we can run a generic collocation code by specifying the scheme at the end of the command. For example, to use a trapezoid scheme:
we can enter:
orb_transfer 3.986004418e5 1500 6000 20 6678 0.5 5 500 2 [0,1] [[0,0],[0.5,0.5]] [0.5,0.5] traj_plot
Implementation
My implementation of this boundary value problem solver can be found here: C++ Code
The resulting plots can be displayed with:
Semantic Segmentation¶
Coming Soon
Semantic segmentation instructions will be added soon.
Uploading an Image¶
The first step is to upload an image to be analyzed. To do that expand the menu at the top right corner, above the console screen and upload an image. Give a name to this image, for example image1. This will be the handle for working with the image from the console.
After uploading the file we can display the uploaded image with gdisplay image1.