In order to successfully complete this assignment, you must follow all the instructions in this notebook and upload your edited ipynb file to D2L with your answers on or before 11:59pm on Wednesday November 25th.
BIG HINT: Read the entire homework before starting.
Homework 5: Working with Linear Systems¶
In this homework we are going to explore many different problem examples. These examples are in the style of your bi-weekly quizzes and are provided as a way to practice.
Note, this is also a similar format to your final exam in this class.
Question 1 - (25 points) Cryptography¶
In this problem we are going to write some code to encrypt and decrypt message using linear algebra and inverse functions. We will represent the message as a matrix $m$ and the encrypted (secret) message as $s$. We will also define an encryption matrix $A$ and a decryption matrix $B$ such that:
$$s=Am$$$$m=Bs$$Step 1 (letters to numbers): The first step is to convert our message string into a list of numbers and be able to take a list of numbers and convert it back into a string. For example; A=0, B=1, C=2, etc. The following functions use an existing letter to number encoding called ASCII so it includes uppercase letters, lowercase letters, numbers and many symbols.
# Here are some libraries you may need to use
%matplotlib inline
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
sym.init_printing(use_unicode=True)
def str2numlist(s):
"""Convert a string to a list of numbers"""
n = [ord(c)-ord('A') for c in s]
return n
def numlist2str(n):
"""Convert list of numbers to a string"""
s = [chr(c+ord('A')) for c in n]
s = ''.join(s)
return s
The following code demonstrates how these functions are used:
message='abcdeABCDE'
num_message = str2numlist(message)
print(num_message)
message = numlist2str(num_message)
print(message)
Step 2 (message as a matrix): Now we need our number list to be in a matrix format. The size of our message matrix depends on the size of our square encryption matrix $A$. For this problem lets assume $A$ is a $4x4$ square matrix. This means that our message matrix should be $4 x l$, where length $l$ changes in size depending on the size of the message.
The following functions use the str2numlist and numlist2str functions from above to convert our message strings into a message matrix $m$ and our message matrix $m$ back to a string.
def message2mat(message, spacing=4):
"""Add extra spaces to message and turn into numpy matrix of size 4x(N%4)"""
message_len = len(message)
extra_spaces = spacing - message_len % spacing
message = message.ljust(len(message) + extra_spaces)
num_message = str2numlist(message)
message_mat=np.matrix(np.reshape(num_message, (spacing,-1)))
return message_mat
def mat2message(message_mat):
try:
"""Unroll a numpy matrix into a list of numbers"""
# found this here https://stackoverflow.com/questions/5183533/how-to-make-list-from-numpy-matrix-in-python
num_message = np.round(np.array(message_mat).reshape(-1,)).astype(int).tolist()
message = numlist2str(num_message)
message = message.strip()
except:
raise Exception('Invalid Message Matrix!!!!')
return message
The following code demonstrates how these functions are used:
m = message2mat(message)
m
message = mat2message(m)
message
Step 3 (encrypt message): Now that we have a message as a matrix we can multiply that matrix by our encryption matrix $A$ which is defined as follows:
A = np.matrix([[ -18, 9, 10, 11], [ 7, -4, 3, 2], [14, -4, 6, 6], [9, 80, 3, -9]])
sym.Matrix(A)
**Question 1.a:** Left Multiply $A$ to the message matrix $m$ and store the results in a new matrix $s$.
#put your code here
from answercheck import checkanswer
checkanswer.detailedwarnings = False
checkanswer.matrix(s,"f5731920bdc98a9255a5cecc3575f331");
Step 4 (Decrypt the message): In order to decrypt your message your friends will already need the decryption matrix $B$.
**Question 1.b:** Generate the decryption matrix $B$
#Put your code here
# YOUR CODE HERE
raise NotImplementedError()
**Question 1.c:** Now left multiply the Decryption matrix $B$ to our encrypted message $s$ and pass the resulting matrix (decripted_message_mat) into mat2message to see the original message
#Put your code here
decripted_message = mat2message(decripted_message_mat)
decripted_message
message == decripted_message
**Question 1.d:** The following secret_message was encrypted using matrix $A$ from above. What is the message?
secret_message = [[ 1124, -202, 361, -389, 1912, 527, -283, 1975, -456, -1068, -1232, -1287, 741],
[ 178, 122, 331, 238, -193, 352, 712, -213, 291, 514, 30, 50, -479],
[ 648, 480, 868, 594, -120, 948, 1360, -138, 702, 830, 166, 162, -824],
[ 3910, 2718, 2882, 4565, 3048, 2800, -2334, 3543, 4619, -1803, 4109, 3199, 4018]]
# YOUR CODE HERE
raise NotImplementedError()
**Question 1.e:** Compose your own secret message (minimum 10 characters long) to the instructor using the above encoding matrix $A$. Make sure you print your original message string, the original message matrix, and the secret message matrix. Include code to check and verify that the message can be decoded using the $B$ matrix.
# YOUR CODE HERE
raise NotImplementedError()
This is a common toy robot arm. It has one motor at the base that rotates the arm around the z-axis and three arm joints that allow you to position the gripper by rotating around the y-axis. Let's model our robot using the following system diagram (note the robot is in its "zero" position):
NOTE: The origin for this robot is located at the base of the robot where it attaches to the ground and is in-line with the first, second and third joints. The x-direction goes from the origin to the right and the z-axis goes from the origin upwards.
The following code uses sympy to represent the $J_5$ transformation, which is a coordinate change from the gripper (end effector) to the fourth joint (moves with the forth joint), and the $J_4$ transformation that goes from the fourth joint to the third joint and takes into account the rotation of the fourth joint and moves with the third joint.
# Here are some libraries you may need to use
%matplotlib inline
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
sym.init_printing(use_unicode=True)
q1,q2,q3,q4,v1,v2,h1,h2 = sym.symbols('\Theta_1, \Theta_2, \Theta_3, \Theta_4, V_1, V_2,H_1,H_2', negative=False)
J5 = sym.Matrix([[1, 0, 0, h2], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
J5
J4 = sym.Matrix([[sym.cos(q4), 0, sym.sin(q4), h1], [0, 1, 0, 0], [-sym.sin(q4), 0, sym.cos(q4), 0], [0, 0, 0, 1]])
J4
✅ **Question 2.a:** Construct two ($4 \times 4$) transformation matrices called $J_3$, and $J_2$. Each transform represents a coordinate change in the toy robot. $J_3$ is the coordinate change from the third joint to the second joint (this takes into consideration the rotation around the third joint, denoted $\theta_3$ degrees); $J_2$ is the coordinate change from the second joint to the first joint (this takes into consideration the rotation around the second joint, denoted $\theta_2$ degrees). Represent your matrices using sympy and the provided symbols:
# YOUR CODE HERE
raise NotImplementedError()
✅ **Question 2.b:** Construct a joint transformation matrix called $J_1$ that transforms from the joint 2 coordinate system to the world coordinate system. This transformation takes into consideration the rotation around the first joint by $\theta_1$ degrees. Represent your matrix using sympy and the provided symbols:
# YOUR CODE HERE
raise NotImplementedError()
✅ **Question 2.c:** The following code is modified from the robotics in-class assignment. Rewrite the joint transformation matrices ($J1, J2$ and $J3$) from questions 2.a and 2.b as numpy matrices with discrete values (instead of symbolic). Replace the identity matrices below with your transformations and use this to simulate the robot and verify your answers:
# ROBOT SIMULATOR
%matplotlib inline
import numpy as np
import math
import matplotlib.pylab as plt
from ipywidgets import interact
from mpl_toolkits.mplot3d import Axes3D
def Robot_Simulator(a1=0,a2=-0,a3=0, a4=0, show=True):
#Convert input joint angles to radians
theta1 = a1/180 * np.pi
theta2 = a2/180 * np.pi
theta3 = a3/180 * np.pi
theta4 = a4/180 * np.pi
H1 = 4
H2 = 2
V1 = 3
V2 = 4
#Define your transfomraiton matrices here.
J1 = np.matrix([[1, 0, 0, 0 ],
[0, 1, 0, 0 ],
[0, 0, 1, 0 ],
[0, 0, 0, 1 ]])
J2 = np.matrix([[1, 0, 0, 0 ],
[0, 1, 0, 0 ],
[0, 0, 1, 0 ],
[0, 0, 0, 1 ]])
J3 = np.matrix([[1, 0, 0, 0 ],
[0, 1, 0, 0 ],
[0, 0, 1, 0 ],
[0, 0, 0, 1 ]])
J4 = np.matrix([[np.cos(theta4),0,np.sin(theta4),H1],
[0, 1, 0, 0 ],
[-np.sin(theta4),0,np.cos(theta4),0],
[0, 0, 0, 1]])
J5 = np.matrix([[1, 0, 0, H2 ],
[0, 1, 0, 0 ],
[0, 0, 1, 0],
[0, 0, 0, 1]])
#Make the rigid end effector
p = np.matrix([[0.5,-0.5,0, 1],
[0,-0.5,0,1],
[0,0.5,0,1],
[0.5,0.5,0,1],
[0,0.5,0,1],
[0,0,0,1]]).T
#Propogate points though the simulation
p = np.concatenate((J5*p, np.matrix([0,0,0,1]).T), axis=1 )
p = np.concatenate((J4*p, np.matrix([0,0,0,1]).T), axis=1 )
p = np.concatenate((J3*p, np.matrix([0,0,0,1]).T), axis=1 )
p = np.concatenate((J2*p, np.matrix([0,0,0,1]).T), axis=1 )
p = np.concatenate((J1*p, np.matrix([0,0,0,1]).T), axis=1 )
if show:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(p[0,:].tolist()[0],(p[1,:]).tolist()[0], (p[2,:]).tolist()[0], s=20, edgecolors='r')
ax.scatter(0,0,0, s=20, facecolors='r', edgecolors='r')
ax.scatter(5,0,0, s=200, facecolors='r', edgecolors='r')
ax.plot(p[0,:].tolist()[0],(p[1,:]).tolist()[0], (p[2,:]).tolist()[0])
ax.set_xlim([-6.5,6.5])
ax.set_ylim([-6.5,6.5])
ax.set_zlim([0,10.5])
ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')
ax.set_zlabel('z-axis')
plt.show()
return p[:,0].T
target = interact(Robot_Simulator, a1=(-180,180), a2=(-180,180), a3=(-180,180), a4=(-180,180));
✅ **Question 2.d:**: Given a point ($P_g$) in gripper coordinates and the above transformation matrices ($J_1,J_2,J_3,J_4,J_5$). Write down the equation to transform $P_g$ into world coordinates $P_w$.
YOUR ANSWER HERE
✅ **Question 2.e:**: Now write the equation to transform $P_w$ into gripper coordinates.
YOUR ANSWER HERE
Question 3 - (25 points) Oscillating Systems¶
The motion of a weight attached to a spring is governed by Hook's Law:
$$ \text{tension} = k \times \text{extension}$$where $k$ is a constant of the spring.
Consider the system shown in the next figure.
$x_1$, $x_2$ and $x_3$ are the displacements of the three balls shown. Each ball sits on a perfectly smooth surface (i.e. no friction) and has equal mass ($m$). The application of Hooke's Law on each ball gives the following second order system of equations:
$$m \ddot x_1 = -k_1x_1 + k_2(x_2 - x_1) $$
$$m \ddot x_2 = -k_2(x_2 - x_1)+ k_3(x_3 - x_2)$$
$$m \ddot x_3 = - k_3(x_3 - x_2) + k_4x_3$$
Given $k_1 = 1$, $k_2 = 2$, $k_3 = 3$, $k_4 = 3$ and $m=2.3$ answer the following questions:
# Here are some libraries you may need to use
%matplotlib inline
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
sym.init_printing(use_unicode=True)
**Question 3.a:** The above system can be rewritten in matrix form mat such that:
$$ \ddot{x} = A x$$$$\text{or}$$$$ \left[\begin{matrix} \ddot{x_1} \\ \ddot{x_2} \\ \ddot{x_3} \end{matrix} \right] = \left[ \begin{matrix} -\frac{k_1 + k_2}{m} & \frac{k_2}{m} & 0 \\ \frac{k_2}{m} & \frac{2k3}{m} & \frac{k_3}{m} \\ 0 & \frac{k_3}{m} & \frac{k_4-k_3}{m} \end{matrix} \right] \left[\begin{matrix} x_1 \\ x_2 \\ x_3 \end{matrix} \right] $$Substitute the above variables and write $A$ as a numpy matrix:
##Put your answer here
from answercheck import checkanswer
checkanswer.detailedwarnings = False
checkanswer.matrix(A,"7c04bd6318542f6ea78912c9a062ea9e");
**Question 3.b:** The above matrix is symmetric and therefore diagonalizable. Calculate $D$ and $C$ such that $D$ is diagonal and $A = CDC^{-1}$:
# YOUR CODE HERE
raise NotImplementedError()
**Question 3.c:** Show that indeed $A=CDC^{-1}$ by comparing the original matrix $A$ to the matrix resulting from multiply $C$ and $D$ from the question 3.b.
# YOUR CODE HERE
raise NotImplementedError()
**Question 3.d:** Now we can rewrite the above equation in terms of $C$ and $D$: $$ \ddot{x} = CDC^{-1} x$$ Which can be transformed into: $$ C^{-1}\ddot{x} = DC^{-1} x$$ Lets introduce a new coordinate system defined by: $$ x' = C^{-1}x$$ $$\text{and}$$ $$ \ddot{x}' = C^{-1}\ddot{x}$$ Rewriting the above equations in the new coordinate system we get:
$$ \ddot{x}' = D x'$$Explain why the system of differential equations in the new coordinate system easier to solve than our previous one? i.e. what is the purpose of the coordinate system change?
YOUR ANSWER HERE
**Question 3.e:** The solutions ($x'$) to the above second order system is of the form:
$$x_n'= b_n \cos \left( \sqrt{d_{nn}} t + \gamma_n \right)$$Where $n$ is our variable index, $d_{nn}$ is the corresponding diagonal component of the $D$ matrix, $t$ is time, and $b_n, \gamma_n$ are constants scalars in $R$.
Explain how to you calculate the final solution ($x$) in the original coordinate system? Note, you do not need to do the calculation. Just explain how the final solution is calculated.
YOUR ANSWER HERE
Question 4 - (25 points) Forward Kinematics #2¶
Consider the crane depicted in the following image.
This style of crane is often called a "crawler" because it has tank like treads to "crawl" around. For this question we will ignore the treads and focus on its actuators. It has; one motor that rotates the cabin around the z-axis; two "arms" which moves up-and-down; and two "Hoist" lines, one connected to the end of the first "arm" (which is also called the boom) and one connected to the second "arm" (which is also called the jib). Each of the hoist lines move up and down due to gravity and are therefore always vertical.
This system can be written as the following diagram:
The origin for this system is located at the base of treads and is in-line with the first joint. The x-direction goes from the origin to the right and the z-axis goes from the origin upwards.
The following code simulates the forward kinematics of the system.
%matplotlib inline
import matplotlib.pylab as plt
from ipywidgets import interact
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import copy
def Robot_Simulator(theta1=0,theta2=45,theta3=-10, D1=1, D2=1):
#Convert from degrees to radians
q1 = theta1/180 * np.pi
q2 = theta2/180 * np.pi
q3 = theta3/180 * np.pi
V1 = 2
H1 = 1
A1 = 4
A2 = 2
J1 = np.matrix([[np.cos(q1), -np.sin(q1), 0, 0 ],
[np.sin(q1), np.cos(q1), 0, 0 ],
[0, 0, 1, V1],
[0, 0, 0, 1]])
J2 = np.matrix([[np.cos(q2), 0, -np.sin(q2), H1 ],
[0, 1, 0, 0],
[np.sin(q2), 0, np.cos(q2), 0 ],
[0, 0, 0, 1]])
J3 = np.matrix([[1, 0, 0, A1 ],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
J4a = np.matrix([[np.cos(-q2), 0, -np.sin(-q2), 0 ],
[0, 1, 0, 0],
[np.sin(-q2), 0, np.cos(-q2), 0 ],
[0, 0, 0, 1]])
J4b = np.matrix([[1, 0, 0, 0 ],
[0, 1, 0, 0],
[0, 0, 1, -D1 ],
[0, 0, 0, 1]])
J5 = np.matrix([[np.cos(q3), 0, -np.sin(q3), 0 ],
[0, 1, 0, 0],
[np.sin(q3), 0, np.cos(q3), 0 ],
[0, 0, 0, 1]])
J6 = np.matrix([[1, 0, 0, A2 ],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
J7 = np.matrix([[np.cos(-q3-q2), 0, -np.sin(-q3-q2), 0 ],
[0, 1, 0, 0],
[np.sin(-q3-q2), 0, np.cos(-q3-q2), 0 ],
[0, 0, 0, 1]])
J8 = np.matrix([[1, 0, 0, 0 ],
[0, 1, 0, 0],
[0, 0, 1, -D2 ],
[0, 0, 0, 1]])
p = np.matrix([[-0.25,0,-0.25, 1], [-0.25,0,0,1], [0.25,0,0, 1], [0.25,0,-0.25,1],[0.25,0,0, 1], [0,0,0,1]]).T
p2 = copy.deepcopy(p)
#Propogate hoist points though the simulation
p = np.concatenate((J8*p, np.matrix([0.0,0,0,1]).T), axis=1 )
p = np.concatenate((J7*p, np.matrix([0.0,0,0,1]).T), axis=1 )
p = np.concatenate((J6*p, np.matrix([0.0,0,0,1]).T), axis=1 )
p = np.concatenate((J5*p, np.matrix([0,0,0,1]).T), axis=1 )
p = np.concatenate((J3*p, np.matrix([0,0,0,1]).T), axis=1 )
p = np.concatenate((J2*p, np.matrix([0,0,0,1]).T), axis=1 )
p = np.concatenate((J1*p, np.matrix([0,0,0,1]).T), axis=1 )
p2 = np.concatenate((J4b*p2, np.matrix([0,0,0,1]).T), axis=1 )
p2 = np.concatenate((J4a*p2, np.matrix([0,0,0,1]).T), axis=1 )
p2 = J3*p2
p2 = J2*p2
p2 = J1*p2
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(p[0,:].tolist()[0],(p[1,:]).tolist()[0], (p[2,:]).tolist()[0], s=20, edgecolors='r')
ax.scatter(p2[0,:].tolist()[0],(p2[1,:]).tolist()[0], (p2[2,:]).tolist()[0], s=20, edgecolors='r')
ax.scatter(0,0,0, s=20, facecolors='r', edgecolors='r')
ax.plot(p[0,:].tolist()[0],(p[1,:]).tolist()[0], (p[2,:]).tolist()[0], color='blue')
ax.plot(p2[0,:].tolist()[0],(p2[1,:]).tolist()[0], (p2[2,:]).tolist()[0], color='blue')
ax.set_xlim([-5,5])
ax.set_ylim([-5,5])
ax.set_zlim([0,6])
ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')
ax.set_zlabel('z-axis')
plt.show()
target = interact(Robot_Simulator, theta1=(-180,180), theta2=(-45,160), theta3=(-45,160), D1=(0.25,4.0), D2=(0.25,4.0)); ##TODO: Modify this line of code
**Question 4.a:** Matrix J2 represents the angle (q2) of rotation of the first boom and matrix J5 represents the angle (q3) of rotation of the jib. Describe the transformation matrix J7 in the above code. Explain why is it using -q3-q2 as an angle for it's rotation?
YOUR ANSWER HERE
**Question 4.b:** Define a new transformation matrix (T1) in terms of J1-J8 that will convert from hoist coordinates at the end of D2 down to the world coordinates at the base of the crane?
YOUR ANSWER HERE
**Question 4.c:** Define another new transformation matrix (T2) in terms of J1-J8 that will convert from world coordinates at the base of the crane to the end of D1.
YOUR ANSWER HERE
**Question 4.d:** Define a third transformation matrix T3 in terms of J1-J8, T1, T2 that will convert from coordinates at the end of D1 to the end of D2.
YOUR ANSWER HERE
**Question 4.e:** The code at line 64 is defining a matrix $P$. What does this matrix represent?
YOUR ANSWER HERE
Congratulations, we're done!¶
Turn in your assignment using D2L no later than 11:59pm on the day of class. See links at the end of this document for access to the class timeline for your section.
Written by Dirk Colbry and Ming Yan, Michigan State University
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.