Projects

Below are many of the projects that I have worked on. Those whose code is not proprietary are accompanied with a link to their respective repos.

Table of Contents

Project Description Repo
Star Tracker Testing for a low-cost cubesat star-tracker
Rigid Body Dynamics Simulation Rigid body dynamics and control system simulation for the Cal Poly Cubesat Lab
Stickcube Simulation, building and testing of self-balancing pendulum for ADCS testing purposes Repository
EIC Crab Cavity LLRF Simulations Simulation and design specifications for the EIC LLRF
Cyclist Energy Optimization, MCM 2022 Optimation with SQP of road cyclist bioenergy for winning entry in 2022 MCM Repository
Fermilab Muon Beamline Analysis Simulation and data analysis of Fermilab's Muon Campus M4 beamline Repository
Mycology Data Analysis Clustering and topological data analysis on Cal Poly Mycology club data Repository

Cal Poly Cubesat Laboratory - Star Tracker

Star trackers are a powerful attitude determination device that is used in many commercial satellite missions. Until recently, star trackers have been mostly out of reach of low-cost cubesat missions due to their cost. Using software and hardware compatible with CPCL missions, we were able to create a star tracker achieving 89 percent accuracy in star identification for $250, compared with commercial star trackers nearing $30,000. The Blue Canyon Standard NST (pictured on right) is capable of 6 arcsec attitude knowledge (much higher fidelity than required for a low-budget mission) but is far beyond the monetary reach of a typical low-budget mission.

The project involved integration and testing of star tracker software written by Luca Merlo Paula Soares. His star tracker software was written in C for use in the CPCL modified linux kernel.

My part in the project was to assist in testing the star tracker's and verifying the attitude algorithms implemented. This involved writing python code to visualize and compare quaternion outputs from the software as well as checking the code's implementation of QUEST (the chosen algorithm to generate a quaternion). Quaternions are 4-dimensional numbers that represent a 3-dimensional rotation so visualizing them has to be done carefully. To do this, I utilized mplot3d (part of matplotlib) to visualize 3d vectors rotated by the quaternion, as well as computing quaternions from the same data using different algorithms, and comparing the hyperspherical distance between them. This must be done carefully as there are infinitely many paths along a hypersphere between 2 quaternions. The TRIAD algorithm was used to generate test quaternions from the same data. I also worked to verify the spherical projection being done in the star tracker software via a 3-dimensional visualization of the vectors and projected camera image.

Hyperspherical distance was computed via the following

def hyp_arc_dist(q1,q2):
	# computes the hyperspherical arc-length (great circle distance) between quaternion 1 and 2 (given by q1 and q2)
	# quaternions are assumed to be given in the (real,i,j,k) convension

	# compute cosine of angle between then
	cospsi = q1[0]*q2[0]+q1[1]*q2[1]+q1[2]*q2[2]+q1[3]*q2[3]

	# obtain angle
	psi = math.acos(cospsi)

	# dont need to scale by the radius because these are normalize (radius 1)
	return 1*psi

This work lead to a published paper in the the American Institute of Aeronautics and Astronautics (AIAA) Journal, which can be found here.

Cal Poly Cubesat Laboratory - Large Scale Rigid Body Dynamics Simulation for ADCS

The attitude determination and control system (ADCS) is historically a difficult system to design and implement. While, the Cal Poly Cubesat Laboratory has historically used simplified, passive systems, increasing mission complexity has called for a more robust design framework for the CPCL ADCS.

This project involved working off of the MATLAB / Simulink framework developed by Liam Bruno as part of his masters thesis in aerospace engineering, which can be found here. The code features an adjustable orbit propagator via the equation \[\mathbf{a}=\frac{\mu}{|\mathbf{R}|^3}\mathbf{R}\] a simulation of the rigid-body dynamics of cubesat via the Euler equation \[\mathbf{I}\mathbf{\omega} + \mathbf{\omega}\times(\mathbf{I}\mathbf{\omega})=\mathbf{M}\] a P.I.D. controller and tuner as well as several control modes.

The goal of this project was 2-fold: to add additional control modes in order to accommodate upcoming missions and to add a simulation of the cubesat power system to the simulink. In many missions, the solar cells are not distributed evenly around the satellite so the power system depends heavily on the attitude of the satellite.

For this project I implemented a control mode for the \(\dot{B}\) algorithm, a detumbling algorithm involving only magnetometers and magnetorquers. I also worked to implement an active magnetorquer-only control algorithm. Additionally, I added power system simulation capability. The simulation now outputs power system information through the different control modes. The tool was used for the CPCL mission AMDROHP to compare different power budgets under different control modes.

Cal Poly Cubesat Laboratory / ETOILES - Stickcube

The attitude determination and control system (ADCS) is often the most difficult to test because the microgravity conditions in space are difficult to simulate on earth without expensive equipment. In an attempt to remedy this for CPCL, Stickcube was started to be a testing apparatus for attitude determination hardware and software as well as control algorithms. The apparatus would involve an inverted pendulum, controlled with two reaction wheels. The system was designed to use an arduino uno mounted to the top of the pendulum with the BNO055 IMU unit and 2 Pololu 9.7:1 metal gearmotors. A P.I.D. controller was programmed onto the arduino along with a complementary filter to obtain the IMU's absolute orientation.

I was the lead for the project over its nearly 3 year lifetime. For the project, several cad models were made and verified for manufacturability. The wheels were manufactured out of aluminum 6061. To get the wheel and controller specifications, a full simulation was created and tested in simulink. The simulation included the rigid body dynamics of the pendulum, the motor and reaction wheel dynamics and the P.I.D. controller. Additional simulation work was done in a python simulation to better account for the pendulum's angular acceleration in the attitude readings.

The complementary filter was discretized using the bilinear transform from a continuous-time s-domain transfer function to a discrete time z-domain transfer function. The following difference equation was then obtained for its implementation in the code \[\Theta_{n} = \frac{1}{2+\omega_{c}T}(T*(u_{n}+u_{n-1})+\Theta_{n-1}*(\omega_c T - 2)) \] where \(\omega_c\) is the complementary filter cutoff frequency, \(T\) is the sample time, and \(u_n=g_n+\omega_c a_n\) where \(g_n\) is the IMU's current gyroscope reading and \(a_n\) is the IMU's current accelerometer reading. Similarly the bilinear transform was used to convert the continuous time P.I.D. controller equation \[ C(t) = k_p*e(t) + k_d*e'(t) + k_i*\int_{0}^{t}e(\tau)d\tau \] to the difference equation \[ c_n = \frac{T}{2}(a*e_n + b*e_{n-1} + c*e_{n-2} + c_{n-2}) \] where \[ a=\frac{2k_p}{T}+\frac{4k_d}{T^2}+k_i, b=2k_i - \frac{8k_d}{T^2}, c=\frac{4k_d}{T^2}+k_i-\frac{2k_p}{T} \] and \(e_n\) is the error in the angle outputted by the complementary filter.

These algorithms were implemented in the code architecture of a finite state machine, written in C++ for arduino. To perform tests on the filter parameters, a test script was written in python to read outputs of the filter and visualize the angles output with matplotlib. This test script was used to determine the optimal cutoff frequency, \(\omega_c\). Work was also done to design a desaturation manuever to desaturate the system's reaction wheels. An additional servo motor apparatus was designed to be placed on top of the pendulum, which would adjust the servo according to the law \[\theta = \text{asin}\left(\frac{K\omega}{mgl}\right) \] where \(K\) is the gain of the desat law, \(\omega\) is the wheel's current speed, \(m\) is the mass of servo arm, \(l\) is the length of the arm and \(g\) is the acceleration due to gravity. Additional work was also done to implement a Kalman filter to replace the complementary filter, taking in the gyroscope and accelerometer data.

This work was originally done in the Cal Poly Cubesat Laboratory, but was later moved to the cubesat lab ETOILES under the guidance of Dr. Pauline Faure. Due to staff changes and other funding reasons, we were not able to finish the project. But the code still remains and can be found at the following repo.

Repository

Cal Poly Physics Department Under Dr. Themis Mastoridis - EIC Crab Cavity LLRF Simulations and Specifications

In a circular particle accelerator, crab cavities are electromagnetic resonators with a resonance in the RF spectrum. They are used to tilt the beam before and after the point of interaction in order to increase beam overlap and thus increase beam luminosity. The crabbed beam interaction can be seen visually via the above picture (credit: Derong Xu, Study of Harmonic Crab Cavity in EIC Beam-Beam Simulations, 2022). The Electron Ion Collider (EIC) will utilize crab cavities for a 25 mrad tilt before and after the interaction point, increasing the luminosity by an order of magnitude.

Transient Beam Loading

The control system for these cavities must be able to react to the "beam loading", the field disturbances caused by charged particles coming through the cavity every turn. To test the different controller architectures and parameters, a simulation was developed in simulink including the cavity-beam interactions and the Low-level RF (LLRF) controller. The performance of different controller architectures were judged based off transient cavity voltage, beam position and controller power usage.

Myself and two other Cal Poly students, Trevor Hidalgo and Matti Toivola, developed the simulink simulation based on a previous accelerator cavity simulink developed by Dr. Mastoridis. First studied was the RF controller, a P.I. controller, which was tuned based on the above metrics.

Also studied was the one-turn feedback controller. This controller samples the beam current and acts with a comb-like frequency response, exhibiting high gain at the betatron sidebands of the revolution harmonics.

Additionally studied was the global controller, a low-bandwidth controller that samples the sum of the crabbing and uncrabbing cavities and acts on the crabbing. Its intended purpose is for the case of full cavity loss, to lower the other cavity's voltage in order to prevent catastrophe before the beam is dumped.

From this work a Brookhaven National Lab technical note was published, which can be found here. This study is currently being combined with future studies to inform the LLRF controller design.

Impedance Reduction

When the beam passes through the cavity, the wake fields produced by the charged particles interact with future particles further down the beam. Under the right conditions these fields can accumulate and generate unstable modes in the beam's motion. This can lead to the emittance (the area of the beam in phase space) growing, decreasing collisions. The interactions of the beam current with the voltage of the accumulated wake fields is quantified with a frequency-dependent impedance.

Mike Blaskiewicz of Brookhaven National Laboratory has written a simulation in fortran to determine the stability of the beam under these conditions based on a given impedance. Code developed by Dr. Mastoridis in matlab is able to accurately model the cavity's impedance with different controller configurations. My part in this project is to modify Blaskiewicz's fortran code and Dr. Mastoridis' matlab code as well as develop further code in order to integrate the simulations. This involved transferring the complex frequency-dependent impedance data from the matlab code, converting it to a usable format and input it into the fortran code. This work was done mainly using pandas in python.

The integration of these simulations has allowed us to study the stability of the beam-beam interactions under different controller configurations. This will give us further insight into the tradeoff that must be made and the competing restrictions on the controller parameters. This work became the content of my senior project, done with Dr. Mastoridis. That senior project can be found here.

COMAP MCM 2022- Road Cyclist Bioenergy Model and Optimal Control

The COMAP Mathematical Contest in Modeling (MCM) is a yearly 4-day competition in applied mathematics. The competition involves researching and creating a model for a real world system given to you by the competition, analyzing your model and compiling your findings into a 25-page research paper. The 2022 MCM was the third time I participated, and the task was to create a model that can inform a long-distance road cyclist on the optimal strategy for distributing their energy over the course of a race.

To meet this task, my team, consisting of myself, Madison Lytle and Callan Whitney viewed the task as a constrained optimal control problem. The optimization was the functional \[\int_{0}^{L}\frac{1}{v(t)}dt=T_f \] the total time to finish the race (in terms of the rides velocity, \(v(t)\)), subject to the constraint equations \[W'_{cap}-\int_{0}^{T_f}W'_{exp}e^{-\delta(t)(t-u)/\tau_{w}}dt=0, P(t)=F(t)v(t),\] \[ F(t)= C_d(v+v_w)^2 + mg(\sin(\phi)+C_r) + m_e\frac{dv}{dt}, P(t)\leq CP+(P_{max}-CP)W_{exp} \] The constraint equations mainly come from the Skiba model; a dynamic model for human bioenergy. In the model \(W'_{cap}\) and \(W'_{exp}\) represents the anaerobic work capacity and anaerobic work expended, respectively. \(F(t)\) is the force exerted by the rider, in terms of the air resistance (characterized by \(C_d\)), the gravitational force (characterized by the course grade, \(\phi\)) and the rider's mass \(m_e\).

Since this was a highly nonlinear problem, the optimization was intractable in most numerical optimization schemes. The method of sequential quadratic programming for its ability to solve nonlinear constrained optimization problems. This method involves adding a quadratic programming subproblem to a normal Newton's method solver. This method could be easily implemented in matlab using the optimization toolbox. Once implemented, the model was used to optimize over several test courses (pictured above).

The effect of environmental perturbations (wind, turning etc.) was also studied on the model. To generalize the model, as well as reduce the computational cost, we studied courses as concatenations of generic pieces (called course elements) and the optimizations of those course elements. Whether the concatenation of course element optimizations is the same as the optimization of the whole course together was also studied but due to time constraints we were not able to fully determine the conditions in which the optimizations would agree.

Our paper received the designation of finalist in the competition and received the MAA award, one of 24 teams out of about 15,000 to receive an award and one of 3 U.S. teams to receive an award. The paper also led to a presentation at Mathfest 2022, MAA's yearly conference and JMM 2023. Our paper submission can be viewed here. An article was also published on the mathematics department website, which can be found here.

Repository

Fermi National Accelerator Laboratory - Muon Campus Beamline Data Analysis

The Fermi National Accelerator Laboratory (one the DOE national labs) is home to the Muon Campus, a facility featuring a muon beam delivery line, which currently serves the muon G-2 experiment and will eventually serve the Mu2e experiment. A new beamline of the muon campus will be used for this experiment for the first time: the M4 line (pictured below). Comparing the phase-space parameters for the muon beam actually delivered by the M4 line to current models can help inform design changes and optimization. In May of 2022 this beamline was commissioned for the first time.

My project as a S.U.L.I. intern over the summer of 2022 under Dr. Diktys Stratakis was to examine this data from the first commissioning of the beamline and compare it to simulation. The simulation used was G4beamline, a Geant4-based simulation giving the beam phase space at different points along the beamline based on the components of the beamline.

The data was analyzed using python's numpy, pandas, scipy, and GEKKO. Numpy and pandas were used to organize the data and scipy and GEKKO were used for fitting. For each point along the beamline where there was a screen, a Guassian was fit to the intensity over width to give the standard deviation of the particles transverse position in x and y (essentially the beam's width).

Two methods were used to measure the full beam phase-space parameters: the 3-screen method and the quad-scan method. The 3-screen method involved measuring the beam width at 3 different points along the beamline and using the relationship \[ \beta=\beta_0-2d\alpha_0+d^2\gamma_0\] relating the phase-space parameters (the Twiss parameters), \(\beta\), \(\gamma\), \(\alpha\) at the different points along the beamline with the same parameters at the beginning of the 3-screen measured with \(d\), the distance from the beginning of the 3-screen measurements. The quad-scan method used many measurements of the beam width down beam from a single quadrupole magnet whose current is slowly increased. Then using the relationship \[ \sigma^2=A\left(1-\frac{d}{f}\right)^2 - 2Bd\left(1-\frac{d}{f}\right) + Cd^2\] where \(\sigma\) is the beam width, \(f\) is the quadrupole focal length, and \(A,B,C\) can be used to find the Twiss parameters at the quadrupole being modulated.

From these two studies of the data I was able to reconstruct the beam phase space at 5 points along the beamline. Comparing these values to the same phase space parameters from the simulated beamline was indicative of the deviations of the simulation from the beamline. Also from this data I was able to determine the beam emittance, the beam's phase space area, which is a measurement of the spread of the beam. Estimates of the beam emittance over time give us an idea of how significant the nonlinear beam interactions are (as linear beam interactions would predict a constant emittance over time). A paper was written from this work, published in the Journal of Instrumentation which can be found here.

Repository

Mathematical Data Science Final Project - Clustering and TDA for iNaturalist Mushroom Foraging Data

The Cal Poly mycology club has access to data from over 700 muchroom identifications, taken from all over the west coast. All the data is stored with the date, time, location of the mushroom, and species. Many useful pieces of information can be obtained from this data such as the locations of unsually low mushroom sightings.

For the final project of our applied mathematics special topics class in mathematical data analysis, my group consisting of myself, Madison Lytle (from whom the mycology data came from), Conner Fletcher, and Jaxon Greene analyzed this data using two techniques covered in the class: clustering and topological data analysis. Single-linkage hierarchical clustering allowed us to group species together to see how the spacial distribution of the mushrooms aligns with their species distribution. To study the distance scale at which the species distribution most aligns with the species, the minimum linkage threshold, \(r\) was examined in terms of the factor \(\gamma\), derived to measure the difference between the clusters formed and cluster that are random in the mushroom species. The minimum of \(\gamma\) (corresponding to the greatest deviation from randomly distributed species) was found on the order of 10 km. This lead us to conclude that the mushroom location does predict species on the scale of separated ecosystems, but not as much on smaller scales. The topological data analysis allowed us to spot certain holes in the data; places where there was an unusual lack of mushroom sightings.

Repository