{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Building Bicycle-Rider Models with SymBRiM\n", "\n", "Welcome to this tutorial on creating **Sym**bolic **B**icycle-**Ri**der **M**odels using SymBRiM!\n", "\n", "In this notebook, we'll embark on a hands-on journey into SymBRiM, focusing on the practical aspects from a model user's perspective. By the end of this tutorial, you'll have achieved the following learning goals:\n", "\n", "- Create a Carvallo-Whipple bicycle model using SymBRiM.\n", "- Create a custom bicycle model to suit specific needs.\n", "- Extend a bicycle model to include a rider.\n", "- Parameterize a model using the [BicycleParameters](https://github.com/moorepants/BicycleParameters) library.\n", "- Simulate models using a simple `Simulator` utility.\n", "- Visualize a model's behavior with ease using the integrated plotting utility powered by [SymMePlot](https://github.com/TJStienstra/symmeplot).\n", "\n", "Before diving into this tutorial, it's recommended that you familiarize yourself with [the introduction](https://mechmotum.github.io/symbrim/guides/introduction.html) and read through the sections on _Software Overview_ to _BRiM Models_ in [the SymBRiM paper](https://doi.org/10.59490/6504c5a765e8118fc7b106c3). This background information will provide valuable context for what we'll cover here.\n", "\n", "## Tutorial Overview\n", "\n", "This tutorial is structured as follows: We'll begin by working with the default Carvallo-Whipple bicycle model, adhering to [Moore's convention](https://moorepants.github.io/dissertation/eom.html). Our journey will encompass configuring the model in SymBRiM, exporting it, and deriving the equations of motion (EoMs). Subsequently, we'll venture into configuring a new bicycle model, this time incorporating fork suspension. This model will be further extended to include an upper-body rider. Next, we'll parametrize these models, enabling us to conduct multiple simulations. Finally, we'll visualize our results as captivating animations.\n", "\n", "## The Default Bicycle Model\n", "\n", "Our starting point is the Carvallo-Whipple bicycle model, following [Moore's convention](https://moorepants.github.io/dissertation/eom.html). The diagram below illustrates the general configuration of this model:\n", "\n", "
\n", "\n", "For our initial model, we'll build the default Carvallo-Whipple bicycle, which aligns perfectly with the image above:\n", "\n", "
\"Configuration
\n", "\n", "When configuring a bicycle in SymBRiM, each component must be assigned a unique name. Symbol names are derived from component names, so using distinct names is essential. For example, if you would provide `\"wheel\"` as the name for both the front wheel and the rear wheel then the symbol for the radius will also be the same symbol. Additionally, SymBRiM applies conventions for certain components. For instance, rear and front frames follow Moore's convention by default.\n", "\n", "_Note: You can explicitly specify convention usage by employing `RigidRearFrame.from_convention(\"moore\", \"rear_frame\")`._\n", "\n", "**Exercise**: The code below initiates the bicycle configuration. Complete it to match the configuration displayed in the image above.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:38.597586Z", "iopub.status.busy": "2024-10-12T22:55:38.597429Z", "iopub.status.idle": "2024-10-12T22:55:39.580926Z", "shell.execute_reply": "2024-10-12T22:55:39.580348Z" } }, "outputs": [], "source": [ "import warnings\n", "\n", "import sympy as sm\n", "import sympy.physics.mechanics as me\n", "\n", "import symbrim as sb\n", "from symbrim.bicycle import RigidRearFrameMoore, WhippleBicycleMoore\n", "from symbrim.brim import SideLeanSeatSpringDamper\n", "from symbrim.rider import PinElbowTorque, SphericalShoulderTorque" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:39.583027Z", "iopub.status.busy": "2024-10-12T22:55:39.582758Z", "iopub.status.idle": "2024-10-12T22:55:39.587814Z", "shell.execute_reply": "2024-10-12T22:55:39.587306Z" } }, "outputs": [], "source": [ "bicycle_def = sb.WhippleBicycle(\"my_bike\")\n", "assert type(bicycle_def) is WhippleBicycleMoore\n", "bicycle_def.rear_frame = sb.RigidRearFrame.from_convention(\"moore\", \"rear_frame\")\n", "assert type(bicycle_def.rear_frame) is RigidRearFrameMoore\n", "bicycle_def.rear_wheel = sb.KnifeEdgeWheel(\"rear_wheel\")\n", "bicycle_def.rear_tire = sb.NonHolonomicTire(\"rear_tire\")\n", "### BEGIN SOLUTION\n", "bicycle_def.ground = sb.FlatGround(\"ground\")\n", "bicycle_def.front_frame = sb.RigidFrontFrame(\"front_frame\")\n", "bicycle_def.front_wheel = sb.KnifeEdgeWheel(\"front_wheel\")\n", "bicycle_def.front_tire = sb.NonHolonomicTire(\"front_tire\")\n", "### END SOLUTION\n", "assert len(bicycle_def.submodels) == 5\n", "assert len(bicycle_def.connections) == 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the model configured we can now let SymBRiM define the entire model using [define_all](https://mechmotum.github.io/symbrim/_autosummary/symbrim.core.base_classes.ModelBase.html#symbrim.core.base_classes.ModelBase.define_all). This method actually calls five other methods in order, of which the last four were also discussed in the four-bar linkage tutorial:\n", "\n", "1. `define_connections`: All models associate the connections with their required submodels.\n", "2. `define_objects`: Create the objects, such as symbols reference frames, without defining any relationships between them.\n", "3. `define_kinematics`: Establish relationships between the objects' orientations/positions, velocities, and accelerations.\n", "4. `define_loads`: Specifies the forces and torques acting upon the system.\n", "5. `define_constraints`: Computes the holonomic and nonholonomic constraints to which the system is subject.\n", "\n", "A reason why you may want to call each of these methods manually in order is that you would like to intervene. For example, you may want to change a symbol as done below." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:39.589664Z", "iopub.status.busy": "2024-10-12T22:55:39.589380Z", "iopub.status.idle": "2024-10-12T22:55:40.087605Z", "shell.execute_reply": "2024-10-12T22:55:40.087138Z" } }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle rear_{wheel r}$" ], "text/plain": [ "rear_wheel_r" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle front_{wheel r}$" ], "text/plain": [ "front_wheel_r" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle Rr$" ], "text/plain": [ "Rr" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle Rf$" ], "text/plain": [ "Rf" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Manually run all steps\n", "bicycle_def.define_connections()\n", "bicycle_def.define_objects()\n", "display(bicycle_def.rear_wheel.radius, bicycle_def.front_wheel.radius)\n", "# Change the symbol names of the rear and front wheel radius\n", "bicycle_def.rear_wheel.symbols[\"r\"] = sm.Symbol(\"Rr\")\n", "bicycle_def.front_wheel.symbols[\"r\"] = sm.Symbol(\"Rf\")\n", "display(bicycle_def.rear_wheel.radius, bicycle_def.front_wheel.radius)\n", "bicycle_def.define_kinematics()\n", "bicycle_def.define_loads()\n", "bicycle_def.define_constraints()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can export this model to a `System` object from `sympy.physics.mechanics`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:40.090070Z", "iopub.status.busy": "2024-10-12T22:55:40.089614Z", "iopub.status.idle": "2024-10-12T22:55:40.101688Z", "shell.execute_reply": "2024-10-12T22:55:40.101268Z" } }, "outputs": [], "source": [ "system_def = bicycle_def.to_system()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At the moment we have not specified any loads, not even gravity. The simplest way to apply gravity is to use `System.apply_uniform_gravity`. To disturb the bicycle we can apply a lateral force at the saddle, and to control the bicycle a torque actuator at the steer.\n", "\n", "_Note: loads and constraints can almost always be applied after exporting a model to a system object. The only instance when this is not the case is when the kinematics have to be altered in any way, e.g. using auxiliary speeds._" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:40.103566Z", "iopub.status.busy": "2024-10-12T22:55:40.103201Z", "iopub.status.idle": "2024-10-12T22:55:40.112164Z", "shell.execute_reply": "2024-10-12T22:55:40.111608Z" } }, "outputs": [ { "data": { "text/plain": [ "((Force(point=rear_frame_wheel_hub_point, force=0),\n", " Force(point=rear_wheel_masscenter, force=g*rear_wheel_mass*ground_frame.z),\n", " Force(point=rear_frame_steer_hub_point, force=0),\n", " Force(point=front_frame_steer_hub_point, force=0),\n", " Force(point=front_frame_wheel_hub_point, force=0),\n", " Force(point=front_wheel_masscenter, force=front_wheel_mass*g*ground_frame.z),\n", " Force(point=ground_origin, force=g*ground_mass*ground_frame.z),\n", " Force(point=rear_frame_body_masscenter, force=g*rear_frame_body_mass*ground_frame.z),\n", " Force(point=front_frame_body_masscenter, force=front_frame_body_mass*g*ground_frame.z),\n", " Force(point=rear_frame_saddle, force=disturbance(t)*rear_frame_body_frame.y)),\n", " (TorqueActuator(steer_torque(t), axis=rear_frame_body_frame.z, target_frame=rear_frame_body_frame, reaction_frame=front_frame_body_frame),))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "normal = bicycle_def.ground.get_normal(bicycle_def.ground.origin)\n", "g = sm.symbols(\"g\")\n", "disturbance = me.dynamicsymbols(\"disturbance\")\n", "steer_torque = me.dynamicsymbols(\"steer_torque\")\n", "system_def.apply_uniform_gravity(-g * normal)\n", "system_def.add_loads(\n", " me.Force(bicycle_def.rear_frame.saddle.point, disturbance * bicycle_def.rear_frame.wheel_hub.axis)\n", ")\n", "system_def.add_actuators(\n", " me.TorqueActuator(steer_torque, bicycle_def.rear_frame.steer_hub.axis,\n", " bicycle_def.rear_frame.steer_hub.frame, bicycle_def.front_frame.steer_hub.frame)\n", ")\n", "system_def.loads, system_def.actuators" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before forming the EoMs we need to specify which generalized coordinates and speeds are independent and which are dependent." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:40.113966Z", "iopub.status.busy": "2024-10-12T22:55:40.113778Z", "iopub.status.idle": "2024-10-12T22:55:40.122128Z", "shell.execute_reply": "2024-10-12T22:55:40.121728Z" } }, "outputs": [ { "data": { "text/plain": [ "ValueError('\\nThe number of dependent generalized coordinates 0 should be equal to\\nthe number of holonomic constraints 1.\\n\\nThe number of dependent generalized speeds 0 should be equal to the\\nnumber of velocity constraints 5.')" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "try:\n", " system_def.validate_system()\n", "except ValueError as e:\n", " display(e)\n", "system_def.q_ind = [*bicycle_def.q[:4], *bicycle_def.q[5:]]\n", "system_def.q_dep = [bicycle_def.q[4]]\n", "system_def.u_ind = [bicycle_def.u[3], *bicycle_def.u[5:7]]\n", "system_def.u_dep = [*bicycle_def.u[:3], bicycle_def.u[4], bicycle_def.u[7]]\n", "system_def.validate_system()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that everything is define, we can form the EoMs (you can ignore the warnings).\n", "\n", "_Note: The reason to use `\"CRAMER\"` as constraint solver is to prevent zero-divisions that otherwise occur._" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:40.123761Z", "iopub.status.busy": "2024-10-12T22:55:40.123451Z", "iopub.status.idle": "2024-10-12T22:55:43.510470Z", "shell.execute_reply": "2024-10-12T22:55:43.509959Z" } }, "outputs": [], "source": [ "with warnings.catch_warnings():\n", " warnings.simplefilter(\"ignore\")\n", " eoms = system_def.form_eoms(constraint_solver=\"CRAMER\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Congratulations you've just created your first bicycle model using SymBRiM!\n", "\n", "## Modifying a Bicycle\n", "Now that you've got a handle on creating a basic bicycle model, let's start modifying it to include some fork suspension. The best method is to create a new `WhippleBicycle` object to prevent possible clashes. The model we will use is [SuspensionRigidFrontFrame](https://mechmotum.github.io/symbrim/_autosummary/symbrim.bicycle.front_frames.SuspensionRigidFrontFrame.html).\n", "\n", "**Exercise**: Complete the configuration of the bicycle model with fork suspension below." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:43.512487Z", "iopub.status.busy": "2024-10-12T22:55:43.512146Z", "iopub.status.idle": "2024-10-12T22:55:43.516494Z", "shell.execute_reply": "2024-10-12T22:55:43.516109Z" } }, "outputs": [], "source": [ "bicycle = sb.WhippleBicycle(\"bicycle\")\n", "### BEGIN SOLUTION\n", "bicycle.ground = sb.FlatGround(\"ground\")\n", "bicycle.rear_wheel = sb.KnifeEdgeWheel(\"rear_wheel\")\n", "bicycle.rear_frame = sb.RigidRearFrame(\"rear_frame\")\n", "bicycle.front_frame = sb.SuspensionRigidFrontFrame(\"front_frame\")\n", "bicycle.front_wheel = sb.KnifeEdgeWheel(\"front_wheel\")\n", "bicycle.rear_tire = sb.NonHolonomicTire(\"rear_tire\")\n", "bicycle.front_tire = sb.NonHolonomicTire(\"front_tire\")\n", "### END SOLUTION\n", "assert len(bicycle.submodels) == 5\n", "assert len(bicycle.connections) == 2\n", "assert isinstance(bicycle.front_frame, sb.SuspensionRigidFrontFrame)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In theory, you could follow the same process now as before and then you will have formed the EoMs of your (first) custom bicycle model. However, we'll continue by using this bicycle model in a bicycle-rider model.\n", "\n", "## Extending with a Rider\n", "We can also extend the model with an upper-body rider. We will attach a side-leaning upper-body rider model, where the shoulders are modeled as spherical joints and the elbows as pin joints. To make this extension we create instantiate a [bicycle-rider model](https://mechmotum.github.io/brim/_autosummary/brim.brim.bicycle_rider.BicycleRider.html). This instance we'll have to provide both our bicycle model, rider model, and the method to connect them. The image below shows the constitution of a complete rider model (not the one we are modeling currently).\n", "\n", "
\"Configuration
\n", "\n", "Let us start with the rider model.\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:43.518145Z", "iopub.status.busy": "2024-10-12T22:55:43.517843Z", "iopub.status.idle": "2024-10-12T22:55:43.521231Z", "shell.execute_reply": "2024-10-12T22:55:43.520848Z" } }, "outputs": [], "source": [ "rider = sb.Rider(\"rider\")\n", "rider.pelvis = sb.PlanarPelvis(\"pelvis\")\n", "rider.torso = sb.PlanarTorso(\"torso\")\n", "rider.left_arm = sb.PinElbowStickLeftArm(\"left_arm\")\n", "rider.right_arm = sb.PinElbowStickRightArm(\"right_arm\")\n", "rider.sacrum = sb.FixedSacrum(\"sacrum\")\n", "rider.left_shoulder = sb.FlexRotLeftShoulder(\"left_shoulder\")\n", "rider.right_shoulder = sb.FlexRotRightShoulder(\"right_shoulder\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let us combine the bicycle and rider models into a single bicycle-rider model. The image below shows the constitution of a complete bicycle-rider model (not the one we are modeling currently). We will just be adding the two models as submodels and the [FixedSeat](https://mechmotum.github.io/brim/_autosummary/brim.brim.seats.FixedSeat.html) and [HolonomicHandGrips](https://mechmotum.github.io/brim/_autosummary/brim.brim.hand_grips.HolonomicHandGrips.html) connections.\n", "\n", "
\"Configuration
\n", "\n", "**Exercise**: Configure the bicycle-rider model accordingly." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:43.522726Z", "iopub.status.busy": "2024-10-12T22:55:43.522441Z", "iopub.status.idle": "2024-10-12T22:55:43.525734Z", "shell.execute_reply": "2024-10-12T22:55:43.525230Z" } }, "outputs": [], "source": [ "bicycle_rider = sb.BicycleRider(\"bicycle_rider\")\n", "### BEGIN SOLUTION\n", "bicycle_rider.bicycle = bicycle\n", "bicycle_rider.rider = rider\n", "bicycle_rider.seat = sb.FixedSeat(\"seat\")\n", "bicycle_rider.hand_grips = sb.HolonomicHandGrips(\"hand_grips\")\n", "### END SOLUTION\n", "assert len(bicycle_rider.submodels) == 2\n", "assert len(bicycle_rider.connections) == 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To control the bicycle we can make the rider torque-driven by applying torques with the elbows. For this purpose SymBRiM has already dedicated load groups, e.g. [PinElbowTorque](https://mechmotum.github.io/symbrim/_autosummary/symbrim.rider.arms.PinElbowTorque.html).\n", "\n", "**Exercise**: Instantiate a load group for each arm with a unique name and assign them to the variables `left_arm_lg` and `right_arm_lg`. Next, add each load group to the corresponding arm using [add_load_groups](https://mechmotum.github.io/symbrim/_autosummary/symbrim.rider.arms.ArmBase.html#symbrim.rider.arms.ArmBase.add_load_groups)." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:43.527505Z", "iopub.status.busy": "2024-10-12T22:55:43.527198Z", "iopub.status.idle": "2024-10-12T22:55:43.530241Z", "shell.execute_reply": "2024-10-12T22:55:43.529858Z" } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "left_arm_lg = PinElbowTorque(\"left_elbow\")\n", "rider.left_arm.add_load_groups(left_arm_lg)\n", "right_arm_lg = PinElbowTorque(\"right_elbow\")\n", "rider.right_arm.add_load_groups(right_arm_lg)\n", "### END SOLUTION\n", "assert isinstance(left_arm_lg, PinElbowTorque)\n", "assert isinstance(right_arm_lg, PinElbowTorque)\n", "assert len(rider.left_arm.load_groups) == 1\n", "assert len(rider.right_arm.load_groups) == 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the model configured we can let SymBRiM define the model.\n", "\n", "_Note: to simplify the EoMs we can set the yaw and roll angles of the rider with respect to the rear frame to zero._" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:43.531979Z", "iopub.status.busy": "2024-10-12T22:55:43.531616Z", "iopub.status.idle": "2024-10-12T22:55:44.799754Z", "shell.execute_reply": "2024-10-12T22:55:44.799253Z" } }, "outputs": [], "source": [ "bicycle_rider.define_connections()\n", "bicycle_rider.define_objects()\n", "bicycle_rider.seat.symbols[\"yaw\"] = 0\n", "bicycle_rider.seat.symbols[\"roll\"] = 0\n", "bicycle_rider.define_kinematics()\n", "bicycle_rider.define_loads()\n", "bicycle_rider.define_constraints()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Export the model to single `System` instance and assign it to the variable `system_br` and apply a gravitational and a disturbance force like before." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:44.801787Z", "iopub.status.busy": "2024-10-12T22:55:44.801462Z", "iopub.status.idle": "2024-10-12T22:55:44.835562Z", "shell.execute_reply": "2024-10-12T22:55:44.835137Z" } }, "outputs": [ { "data": { "text/plain": [ "((Force(point=front_frame_suspension_stanchions, force=(-front_frame_c*front_frame_u(t) - front_frame_k*front_frame_q(t))*front_frame_body_frame.z),\n", " Force(point=front_frame_suspension_lowers, force=(front_frame_c*front_frame_u(t) + front_frame_k*front_frame_q(t))*front_frame_body_frame.z),\n", " Force(point=rear_frame_body_masscenter, force=g*rear_frame_body_mass*ground_frame.z),\n", " Force(point=rear_frame_saddle, force=0),\n", " Force(point=pelvis_masscenter, force=g*pelvis_mass*ground_frame.z),\n", " Force(point=rear_frame_wheel_hub_point, force=0),\n", " Force(point=rear_wheel_masscenter, force=g*rear_wheel_mass*ground_frame.z),\n", " Force(point=rear_frame_steer_hub_point, force=0),\n", " Force(point=front_frame_steer_hub_point, force=0),\n", " Force(point=front_frame_wheel_point, force=0),\n", " Force(point=front_wheel_masscenter, force=front_wheel_mass*g*ground_frame.z),\n", " Force(point=ground_origin, force=g*ground_mass*ground_frame.z),\n", " Force(point=front_frame_body_masscenter, force=front_frame_body_mass*g*ground_frame.z),\n", " Force(point=torso_masscenter, force=g*torso_mass*ground_frame.z),\n", " Force(point=left_shoulder_int_point, force=0),\n", " Force(point=left_shoulder_int_point, force=0),\n", " Force(point=left_arm_upper_arm_masscenter, force=g*left_arm_upper_arm_mass*ground_frame.z),\n", " Force(point=right_shoulder_int_point, force=0),\n", " Force(point=right_shoulder_int_point, force=0),\n", " Force(point=right_arm_upper_arm_masscenter, force=g*right_arm_upper_arm_mass*ground_frame.z),\n", " Force(point=left_arm_forearm_masscenter, force=g*left_arm_forearm_mass*ground_frame.z),\n", " Force(point=right_arm_forearm_masscenter, force=g*right_arm_forearm_mass*ground_frame.z),\n", " Force(point=rear_frame_saddle, force=disturbance(t)*rear_frame_body_frame.y)),\n", " (TorqueActuator(left_elbow_T(t), axis=left_arm_upper_arm_frame.y, target_frame=left_arm_forearm_frame, reaction_frame=left_arm_upper_arm_frame),\n", " TorqueActuator(right_elbow_T(t), axis=right_arm_upper_arm_frame.y, target_frame=right_arm_forearm_frame, reaction_frame=right_arm_upper_arm_frame)))" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "### BEGIN SOLUTION\n", "system_br = bicycle_rider.to_system()\n", "normal = bicycle_rider.bicycle.ground.get_normal(bicycle_rider.bicycle.ground.origin)\n", "g = sm.symbols(\"g\")\n", "system_br.apply_uniform_gravity(-g * normal)\n", "system_br.add_loads(\n", " me.Force(bicycle_rider.bicycle.rear_frame.saddle.point,\n", " disturbance * bicycle_rider.bicycle.rear_frame.wheel_hub.axis)\n", ")\n", "### END SOLUTION\n", "system_br.loads, system_br.actuators" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Set the independent and dependent generalized coordinates and speeds of `system_br`.\n", "\n", "_Hint: You can use the same division for `bicycle_rider.bicycle.q` and `bicycle_rider.bicycle.u` as before for `bicycle_def`._" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:44.837088Z", "iopub.status.busy": "2024-10-12T22:55:44.836936Z", "iopub.status.idle": "2024-10-12T22:55:44.842669Z", "shell.execute_reply": "2024-10-12T22:55:44.842283Z" } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "b, r = bicycle_rider.bicycle, bicycle_rider.rider\n", "system_br.q_ind = [\n", " *b.q[:4], *b.q[5:], *b.front_frame.q\n", "]\n", "system_br.q_dep = [\n", " b.q[4], *r.left_shoulder.q, *r.right_shoulder.q, *r.left_arm.q, *r.right_arm.q\n", "]\n", "system_br.u_ind = [\n", " b.u[3], *b.u[5:7], *b.front_frame.u\n", "]\n", "system_br.u_dep = [\n", " *b.u[:3], b.u[4], b.u[7], *r.left_shoulder.u, *r.right_shoulder.u, *r.left_arm.u, *r.right_arm.u\n", "]\n", "### END SOLUTION\n", "system_br.validate_system()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Form the equations of motion and use `\"CRAMER\"` as constraint solver." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:55:44.844102Z", "iopub.status.busy": "2024-10-12T22:55:44.843960Z", "iopub.status.idle": "2024-10-12T22:56:07.489647Z", "shell.execute_reply": "2024-10-12T22:56:07.489115Z" } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "with warnings.catch_warnings():\n", " warnings.simplefilter(\"ignore\")\n", " eoms = system_br.form_eoms(constraint_solver=\"CRAMER\")\n", "### END SOLUTION\n", "assert system_br.mass_matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Congratulations! You have just built your (first) bicycle-rider model.\n", "\n", "## Parametrization\n", "While it is possible to manually figure out all parameters and create a dictionary to map all symbols accordingly, SymBRiM integrates support for the [BicycleParameters](https://github.com/moorepants/BicycleParameters) package.\n", "The `BicycleParameters` package is designed to generate, manipulate, and visualize the parameters of the Whipple-Carvallo bicycle model. The [repository](https://github.com/moorepants/BicycleParameters/tree/master/data) provides several measured parametrizations.\n", "For instructions on measuring your own bicycle and/or rider refer to [Moore's dissertation](https://moorepants.github.io/dissertation/physicalparameters.html).\n", "\n", "The code below imports the parametrization library and instantiates a parametrization object, which loads in `\"Jason\"` on the Batavus `\"Browser\"` bicycle.\n", "\n", "_Note: This cell will download some parametrization data if it is not already present._" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:07.491807Z", "iopub.status.busy": "2024-10-12T22:56:07.491467Z", "iopub.status.idle": "2024-10-12T22:56:08.053220Z", "shell.execute_reply": "2024-10-12T22:56:08.052767Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Found the RawData directory: /home/runner/work/symbrim/symbrim/docs/tutorials/data/bicycles/Browser/RawData\n", "Looks like you've already got some parameters for Browser, use forceRawCalc to recalculate.\n", "There is no rider on the bicycle, now adding Jason.\n", "Calculating the human configuration.\n" ] } ], "source": [ "from pathlib import Path\n", "\n", "import bicycleparameters as bp\n", "import numpy as np\n", "from utils import download_parametrization_data\n", "\n", "data_dir = Path.cwd() / \"data\"\n", "if not data_dir.exists():\n", " download_parametrization_data(data_dir)\n", "bike_params = bp.Bicycle(\"Browser\", pathToData=data_dir)\n", "bike_params.add_rider(\"Jason\", reCalc=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get a mapping of symbols to values based on this parametrization we can use [get_param_values](https://mechmotum.github.io/brim/_autosummary/brim.core.base_classes.BrimBase.html#brim.core.base_classes.BrimBase.get_param_values) as shown below." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.054771Z", "iopub.status.busy": "2024-10-12T22:56:08.054620Z", "iopub.status.idle": "2024-10-12T22:56:08.123819Z", "shell.execute_reply": "2024-10-12T22:56:08.123297Z" } }, "outputs": [ { "data": { "text/plain": [ "{rear_frame_body_mass: 9.86,\n", " rear_frame_l_bbx: np.float64(0.40374135119327337),\n", " rear_frame_l_bbz: np.float64(0.22043802153129097),\n", " rear_frame_ixx: 0.6473909751608555,\n", " rear_frame_iyy: 1.3163960124977576,\n", " rear_frame_izz: 0.6390248264511185,\n", " rear_frame_izx: -0.16249625380874022,\n", " rear_frame_d1: 0.963149263487593,\n", " rear_frame_l1: 0.3308144077188117,\n", " rear_frame_l2: -0.07398701275972652,\n", " rear_frame_d4: np.float64(0.4225541286665965),\n", " rear_frame_d5: np.float64(-0.5493321255935064),\n", " front_frame_body_mass: 3.22,\n", " front_frame_ixx: 0.28109686661693817,\n", " front_frame_iyy: 0.245827908036,\n", " front_frame_izz: 0.06782716361036183,\n", " front_frame_izx: 0.00634846730073884,\n", " front_frame_d2: 0.4338396131640938,\n", " front_frame_d3: 0.0705000000001252,\n", " front_frame_l3: -0.07654646159268344,\n", " front_frame_l4: -0.47166687226492093,\n", " front_frame_d6: np.float64(-0.17310563866173195),\n", " front_frame_d7: np.float64(0.29),\n", " front_frame_d8: np.float64(-0.37486632400446973),\n", " rear_wheel_mass: 3.11,\n", " rear_wheel_ixx: 0.0904114316323,\n", " rear_wheel_iyy: 0.152391250767,\n", " Rr: 0.340958858855,\n", " front_wheel_mass: 2.02,\n", " front_wheel_ixx: 0.0883826870796,\n", " front_wheel_iyy: 0.149221207336,\n", " Rf: 0.34352982332,\n", " g: 9.81}" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "constants_def = bicycle_def.get_param_values(bike_params)\n", "constants_def[g] = 9.81 # Don't forget to specify the gravitational constant.\n", "constants_def" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The code below defines a function to quickly retrieve all symbols defined in a model.\n", "Next, we obtain all symbols defined in the bicycle models and if there are still symbols missing in our `constants_def` dictionary." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.125786Z", "iopub.status.busy": "2024-10-12T22:56:08.125318Z", "iopub.status.idle": "2024-10-12T22:56:08.128947Z", "shell.execute_reply": "2024-10-12T22:56:08.128569Z" } }, "outputs": [ { "data": { "text/plain": [ "set()" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "missing_symbols = bicycle_def.get_all_symbols().difference(constants_def.keys())\n", "missing_symbols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Create a `constants_br` dictionary for the `bicycle_rider` model and obtain the missing symbols (`missing_symbols`).\n", "\n", "_Note: you may use the same parametrization `bike_params` as before._" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.130539Z", "iopub.status.busy": "2024-10-12T22:56:08.130240Z", "iopub.status.idle": "2024-10-12T22:56:08.371862Z", "shell.execute_reply": "2024-10-12T22:56:08.371342Z" } }, "outputs": [ { "data": { "text/plain": [ "{front_frame_c,\n", " front_frame_d9,\n", " front_frame_k,\n", " left_elbow_T(t),\n", " right_elbow_T(t),\n", " seat_pitch}" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "### BEGIN SOLUTION\n", "constants_br = bicycle_rider.get_param_values(bike_params)\n", "missing_symbols = bicycle_rider.get_all_symbols().difference(constants_br.keys())\n", "### END SOLUTION\n", "assert len(constants_br) >= 87\n", "missing_symbols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To quickly see what this symbol is used for we can ask the model to get a description for each:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.373491Z", "iopub.status.busy": "2024-10-12T22:56:08.373325Z", "iopub.status.idle": "2024-10-12T22:56:08.794027Z", "shell.execute_reply": "2024-10-12T22:56:08.793506Z" } }, "outputs": [ { "data": { "text/plain": [ "'front_frame_k - Spring stiffness of the front suspension.'" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "'front_frame_d9 - Perpendicular distance from the steer axis to the suspension.'" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "'seat_pitch - Pitch angle of the pelvis w.r.t. the rear frame.'" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "'left_elbow_T(t) - Elbow torque of left_arm'" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "'right_elbow_T(t) - Elbow torque of right_arm'" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "'front_frame_c - Damping coefficient of the front suspension.'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for sym in missing_symbols:\n", " display(f\"{sym} - {bicycle_rider.get_description(sym)}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As these symbols have not been defined by the bicycle parameters data we need to define them manually. However, some of them are inputs, so those we'll specify later on." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.795822Z", "iopub.status.busy": "2024-10-12T22:56:08.795651Z", "iopub.status.idle": "2024-10-12T22:56:08.798705Z", "shell.execute_reply": "2024-10-12T22:56:08.798328Z" } }, "outputs": [], "source": [ "constants_br.update({\n", " bicycle_rider.bicycle.front_frame.symbols[\"d9\"]: constants_br[bicycle_rider.bicycle.front_frame.symbols[\"d3\"]],\n", " bicycle_rider.bicycle.front_frame.symbols[\"k\"]: 19.4E3,\n", " bicycle_rider.bicycle.front_frame.symbols[\"c\"]: 9E3,\n", " bicycle_rider.seat.symbols[\"pitch\"]: -0.7,\n", " g: 9.81 # Just in case you did not specify the gravitational constant yet.\n", "})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Warning**: Above we do forget to take into account the inertia and mass of the legs. A good option would have been computing their contribution and adding that to the rear frame's mass and inertia." ] }, { "cell_type": "markdown", "metadata": { "jupyter": { "outputs_hidden": false } }, "source": [ "## Simulation Bicycle Model\n", "\n", "Now that we have a parametrized model we can simulate it. SymBRiM does not provide a simulation utility (yet), but a separate module has here been provided implementing a simple simulation utility. The code below imports this utility and initializes a simulation object for the bicycle model." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.800280Z", "iopub.status.busy": "2024-10-12T22:56:08.799989Z", "iopub.status.idle": "2024-10-12T22:56:08.813485Z", "shell.execute_reply": "2024-10-12T22:56:08.813099Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import numpy as np\n", "from simulator import Simulator\n", "\n", "sim_def = Simulator(system_def)" ] }, { "cell_type": "markdown", "metadata": { "jupyter": { "outputs_hidden": false } }, "source": [ "Before running the actual simulation we need to specify three things:\n", "- `constants`: A dictionary mapping static symbols to values.\n", "- `initial_conditions`: A dictionary mapping generalized coordinates and speeds to initial values.\n", "- `inputs`: A dictionary mapping input symbols to a function taking the time and state as arguments, i.e. `f(t, x) -> float`.\n", "\n", "The code below sets each of these for the bicycle model, where an initial guess for the rear frame pitch is provided and the rear wheel is given an initial angular velocity. The disturbance inputs is a sine wave that increases in amplitude over time. For the steer input, we'll use a simple proportional controller that applies a torque proportional to the roll rate to stabilize the bicycle." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.815080Z", "iopub.status.busy": "2024-10-12T22:56:08.814768Z", "iopub.status.idle": "2024-10-12T22:56:08.818909Z", "shell.execute_reply": "2024-10-12T22:56:08.818465Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "sim_def.constants = constants_def\n", "sim_def.initial_conditions = {\n", " **{xi: 0.0 for xi in system_def.q.col_join(system_def.u)},\n", " bicycle_def.q[4]: 0.314, # Initial guess rear frame pitch.\n", " bicycle_def.u[5]: -3.0 / constants_def[bicycle_def.rear_wheel.radius], # Rear wheel angular velocity.\n", "}\n", "roll_rate_idx = len(system_def.q) + system_def.u[:].index(bicycle_def.u[3])\n", "max_roll_rate, max_torque = 0.2, 10\n", "sim_def.inputs = {\n", " disturbance: lambda t, x: (30 + 30 * t) * np.sin(t * 2 * np.pi),\n", " steer_torque: lambda t, x: -max_torque * max(-1, min(x[roll_rate_idx] / max_roll_rate, 1)),\n", "}" ] }, { "cell_type": "markdown", "metadata": { "jupyter": { "outputs_hidden": false } }, "source": [ "Now we can initialize the simulator objects, which automatically solves the initial conditions for us.\n", "\n", "_Note: the initialization of a simulator object can take a while, so be patient._" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:08.820286Z", "iopub.status.busy": "2024-10-12T22:56:08.820139Z", "iopub.status.idle": "2024-10-12T22:56:11.622290Z", "shell.execute_reply": "2024-10-12T22:56:11.621787Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "{my_bike_q1(t): 0.0,\n", " my_bike_q2(t): 0.0,\n", " my_bike_q3(t): 0.0,\n", " my_bike_q4(t): 0.0,\n", " my_bike_q6(t): 0.0,\n", " my_bike_q7(t): 0.0,\n", " my_bike_q8(t): 0.0,\n", " my_bike_q5(t): np.float64(0.399680398706867),\n", " my_bike_u4(t): 0.0,\n", " my_bike_u6(t): -8.79871551094032,\n", " my_bike_u7(t): 0.0,\n", " my_bike_u1(t): np.float64(2.9999999999999996),\n", " my_bike_u2(t): np.float64(-1.7142933546584113e-28),\n", " my_bike_u3(t): np.float64(3.057822083862947e-28),\n", " my_bike_u5(t): np.float64(-1.2576551361634348e-16),\n", " my_bike_u8(t): np.float64(-8.732866250175556)}" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sim_def.initialize()\n", "sim_def.initial_conditions" ] }, { "cell_type": "markdown", "metadata": { "jupyter": { "outputs_hidden": false } }, "source": [ "With the simulator initialized we can run the simulation. The code below runs the simulation for 2.5 seconds. A DAE solver is used with a relative tolerance of 1e-3 and an absolute tolerance of 1e-6." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:11.623897Z", "iopub.status.busy": "2024-10-12T22:56:11.623740Z", "iopub.status.idle": "2024-10-12T22:56:16.960851Z", "shell.execute_reply": "2024-10-12T22:56:16.960310Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "np.float64(2.49)" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sim_def.solve(np.arange(0, 2.5, 0.01), solver=\"dae\", rtol=1e-3, atol=1e-6)\n", "sim_def.t[-1]" ] }, { "cell_type": "markdown", "metadata": { "jupyter": { "outputs_hidden": false } }, "source": [ "## Visualizing the Bicycle Simulation Results\n", "\n", "Using the [Plotter](https://mechmotum.github.io/symbrim/_autosummary/symbrim.utilities.plotting.Plotter.html) from `SymBRiM` we can easily visualize the model and create an animation based on the simulation results." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:16.962697Z", "iopub.status.busy": "2024-10-12T22:56:16.962390Z", "iopub.status.idle": "2024-10-12T22:56:20.932568Z", "shell.execute_reply": "2024-10-12T22:56:20.931979Z" }, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "from IPython.display import HTML\n", "from matplotlib.animation import FuncAnimation\n", "from scipy.interpolate import CubicSpline\n", "\n", "from symbrim.utilities.plotting import Plotter\n", "\n", "# Create some functions to interpolate the results.\n", "x_eval = CubicSpline(sim_def.t, sim_def.x.T)\n", "r_eval = CubicSpline(sim_def.t, [[cf(t, x) for cf in sim_def.inputs.values()]\n", " for t, x in zip(sim_def.t, sim_def.x.T)])\n", "p, p_vals = zip(*sim_def.constants.items())\n", "max_disturbance = r_eval(sim_def.t)[:, tuple(sim_def.inputs.keys()).index(disturbance)].max()\n", "\n", "# Plot the initial configuration of the model\n", "fig, ax = plt.subplots(subplot_kw={\"projection\": \"3d\"}, figsize=(8, 8))\n", "plotter = Plotter.from_model(bicycle_def, ax=ax)\n", "plotter.add_vector(disturbance * bicycle_def.rear_frame.wheel_hub.axis / max_disturbance,\n", " bicycle_def.rear_frame.saddle.point, name=\"disturbance\", color=\"r\")\n", "plotter.lambdify_system((system_def.q[:] + system_def.u[:], sim_def.inputs.keys(), p))\n", "plotter.evaluate_system(x_eval(0.0), r_eval(0.0), p_vals)\n", "plotter.plot()\n", "X, Y = np.meshgrid(np.arange(-1, 10, 0.5), np.arange(-1, 3, 0.5))\n", "ax.plot_wireframe(X, Y, np.zeros_like(X), color=\"k\", alpha=0.3, rstride=1, cstride=1)\n", "ax.invert_zaxis()\n", "ax.invert_yaxis()\n", "ax.set_xlim(X.min(), X.max())\n", "ax.set_ylim(Y.min(), Y.max())\n", "ax.view_init(19, 14)\n", "ax.set_aspect(\"equal\")\n", "ax.axis(\"off\")\n", "\n", "fps = 30\n", "ani = plotter.animate(\n", " lambda ti: (x_eval(ti), r_eval(ti), p_vals),\n", " frames=np.arange(0, sim_def.t[-1], 1 / fps),\n", " blit=False)\n", "display(HTML(ani.to_jshtml(fps=fps)))\n", "plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulating the Bicycle-Rider Model\n", "Now that we know how to simulate the bicycle model, we can perform a similar simulation for the bicycle-rider model we have created.\n", "\n", "**Exercise**: Initiate a simulator object for the bicycle-rider model and assign in to `sim_br`." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:20.985791Z", "iopub.status.busy": "2024-10-12T22:56:20.985199Z", "iopub.status.idle": "2024-10-12T22:56:20.988299Z", "shell.execute_reply": "2024-10-12T22:56:20.987791Z" } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "sim_br = Simulator(system_br)\n", "### END SOLUTION\n", "assert sim_br.system is system_br" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Set the constants, initial conditions, and inputs of `sim_br`. For the inputs try using a similar control for the elbow torques as before for the steer torque.\n", "\n", "_Hint: You can access the torques of the load groups from their `symbols` dictionary, e.g. `left_arm_lg.symbols[\"T\"]`._" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:20.990023Z", "iopub.status.busy": "2024-10-12T22:56:20.989739Z", "iopub.status.idle": "2024-10-12T22:56:20.995832Z", "shell.execute_reply": "2024-10-12T22:56:20.995416Z" } }, "outputs": [], "source": [ "### BEGIN SOLUTION\n", "sim_br.constants = constants_br\n", "sim_br.initial_conditions = {\n", " **{xi: 0.0 for xi in system_br.q.col_join(system_br.u)},\n", " bicycle_rider.bicycle.q[4]: 0.314, # Initial guess rear frame pitch.\n", " bicycle_rider.bicycle.u[5]: -3.0 / constants_br[bicycle_rider.bicycle.rear_wheel.radius],\n", " bicycle_rider.rider.left_arm.q[0]: 0.7, # Elbow flexion.\n", " bicycle_rider.rider.right_arm.q[0]: 0.7, # Elbow flexion.\n", " bicycle_rider.rider.left_shoulder.q[0]: 0.5, # Shoulder flexion.\n", " bicycle_rider.rider.right_shoulder.q[0]: 0.5, # Shoulder flexion.\n", " bicycle_rider.rider.left_shoulder.q[1]: -0.6, # Shoulder rotation.\n", " bicycle_rider.rider.left_shoulder.q[1]: -0.6, # Shoulder rotation.\n", " bicycle_rider.bicycle.u[3]: 0.3,\n", "}\n", "roll_idx = system_br.q[:].index(bicycle_rider.bicycle.q[3])\n", "roll_rate_idx = len(system_br.q) + system_br.u[:].index(bicycle_rider.bicycle.u[3])\n", "max_roll_rate, max_torque = 0.2, 10\n", "bound = lambda x: max(-1, min(x, 1))\n", "excitation = lambda x: bound(x[roll_rate_idx] / max_roll_rate)\n", "\n", "sim_br.inputs = {\n", " disturbance: lambda t, x: (30 + 30 * t) * np.sin(t * 2 * np.pi),\n", " left_arm_lg.symbols[\"T\"]: lambda t, x: -max_torque * excitation(x),\n", " right_arm_lg.symbols[\"T\"]: lambda t, x: max_torque * excitation(x),\n", "}\n", "### END SOLUTION\n", "assert len(sim_br.constants) == 92\n", "assert len(sim_br.initial_conditions) == 30\n", "assert len(sim_br.inputs) == 3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Initialize the simulator.\n", "\n", "_Hint: If a runtime error from scipy.fsolve is shown, then it was not able to solve the initial conditions. Recommendable is to update your initial conditions to provide a better guess and use `sim_br.solve_initial_conditions()` without reinitializing the simulator object._\n", "\n", "_Note: Initializing the bicycle-rider simulator takes longer than initializing just the bicycle._" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:56:20.997363Z", "iopub.status.busy": "2024-10-12T22:56:20.997104Z", "iopub.status.idle": "2024-10-12T22:59:35.526433Z", "shell.execute_reply": "2024-10-12T22:59:35.525960Z" } }, "outputs": [ { "data": { "text/plain": [ "array([ 3.00000000e+00, 3.26841150e-33, -5.69472114e-16, 3.00000000e-01,\n", " -8.79871551e+00, -0.00000000e+00, -8.73286625e+00, 0.00000000e+00,\n", " -0.00000000e+00, 1.56914295e-27, -2.22877858e-28, 1.39570538e-27,\n", " -4.37780608e-27, -3.28954997e-28, -1.15750563e-27, -1.75919306e-01,\n", " 2.42600844e+01, 1.40657668e+01, 9.59312566e+00, -5.58384985e+00,\n", " -1.55431223e-15, 7.92697958e-01, -7.88318230e+00, 3.49448718e+01,\n", " 1.44833210e+01, -1.93488172e+01, -1.44833210e+01, 1.93488172e+01,\n", " -7.45538540e+00, 7.45538540e+00])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "### BEGIN SOLUTION\n", "sim_br.initialize()\n", "### END SOLUTION\n", "sim_br.eval_rhs(0, [sim_br.initial_conditions[xi] for xi in system_br.q[:] + system_br.u[:]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise**: Peform a simulation over a time-span of 2.5 seconds using a DAE solver.\n", "\n", "_Hint: If the simulation fails it can be that the controls are unstable causing the bicycle to fall or that something else is wrong. A simple solution is to use an ODE integrator to figure out the controls: `sim_br.solve((0, 2.5), t_eval=np.arange(0, 2.5, 0.01))`. Also, feel free to reduce the disturbance force to make the control easier._" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:59:35.528310Z", "iopub.status.busy": "2024-10-12T22:59:35.527876Z", "iopub.status.idle": "2024-10-12T22:59:37.363534Z", "shell.execute_reply": "2024-10-12T22:59:37.362998Z" } }, "outputs": [ { "data": { "text/plain": [ "np.float64(2.49)" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "### BEGIN SOLUTION\n", "# sim_br.solve(np.arange(0, 2.5, 0.01), solver=\"dae\", rtol=1e-3, atol=1e-6)\n", "sim_br.solve((0, 2.5), t_eval=np.arange(0, 2.5, 0.01))\n", "### END SOLUTION\n", "sim_br.t[-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualizing the Bicycle-Rider Simulation\n", "\n", "The simulation can be visualized using similar code as below.\n", "\n", "**Exercise**: Complete the code below by initializing the [Plotter](https://mechmotum.github.io/brim/_autosummary/brim.utilities.plotting.Plotter.html), optionally adding a vector for the disturbance, and running [lambdify_system](https://mechmotum.github.io/brim/_autosummary/brim.utilities.plotting.Plotter.html#brim.utilities.plotting.Plotter.lambdify_system)." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "execution": { "iopub.execute_input": "2024-10-12T22:59:37.365414Z", "iopub.status.busy": "2024-10-12T22:59:37.365110Z", "iopub.status.idle": "2024-10-12T22:59:42.639572Z", "shell.execute_reply": "2024-10-12T22:59:42.638964Z" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Create some functions to interpolate the results.\n", "x_eval = CubicSpline(sim_br.t, sim_br.x.T)\n", "r_eval = CubicSpline(sim_br.t, [[cf(t, x) for cf in sim_br.inputs.values()]\n", " for t, x in zip(sim_br.t, sim_br.x.T)])\n", "p, p_vals = zip(*sim_br.constants.items())\n", "max_disturbance = r_eval(sim_br.t)[:, tuple(sim_br.inputs.keys()).index(disturbance)].max()\n", "\n", "# Plot the initial configuration of the model\n", "fig, ax = plt.subplots(subplot_kw={\"projection\": \"3d\"}, figsize=(8, 8))\n", "### BEGIN SOLUTION\n", "plotter = Plotter.from_model(bicycle_rider, ax=ax)\n", "plotter.add_vector(disturbance * bicycle_rider.bicycle.rear_frame.wheel_hub.axis / max_disturbance,\n", " bicycle_rider.bicycle.rear_frame.saddle.point, name=\"disturbance\", color=\"r\")\n", "plotter.lambdify_system((system_br.q[:] + system_br.u[:], sim_br.inputs.keys(), p))\n", "### END SOLUTION\n", "plotter.evaluate_system(x_eval(0.0), r_eval(0.0), p_vals)\n", "plotter.plot()\n", "X, Y = np.meshgrid(np.arange(-1, 10, 0.5), np.arange(-1, 4, 0.5))\n", "ax.plot_wireframe(X, Y, np.zeros_like(X), color=\"k\", alpha=0.3, rstride=1, cstride=1)\n", "ax.invert_zaxis()\n", "ax.invert_yaxis()\n", "ax.set_xlim(X.min(), X.max())\n", "ax.set_ylim(Y.min(), Y.max())\n", "ax.view_init(19, 14)\n", "ax.set_aspect(\"equal\")\n", "ax.axis(\"off\")\n", "\n", "fps = 30\n", "ani = plotter.animate(\n", " lambda ti: (x_eval(ti), r_eval(ti), p_vals),\n", " frames=np.arange(0, sim_br.t[-1], 1 / fps),\n", " blit=False)\n", "display(HTML(ani.to_jshtml(fps=fps)))\n", "plt.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Congratulations! You have completed this tutorial." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.10" } }, "nbformat": 4, "nbformat_minor": 4 }