diff --git a/.gitignore b/.gitignore index c4bc088..39bfdb1 100644 --- a/.gitignore +++ b/.gitignore @@ -135,12 +135,309 @@ cython_debug/ # Ignore exercises exercises/ -# LaTeX Garbage -*.aux -*.fdb_latexmk -*.fls -*.log -*.synctex.gz +# Logs +logs/ +log/ +*.xlsx -# LaTeX output -*.pdf \ No newline at end of file +# Drawio backups +*.drawio.bkp +*.drawio.dtmp + +# LaTeX Garbage +## Core latex/pdflatex auxiliary files: +*.aux +*.lof +*.log +*.lot +*.fls +*.out +*.toc +*.fmt +*.fot +*.cb +*.cb2 +.*.lb + +## Intermediate documents: +*.dvi +*.xdv +*-converted-to.* +# these rules might exclude image files for figures etc. +# *.ps +# *.eps +# *.pdf + +## Generated if empty string is given at "Please type another file name for output:" +*.pdf + +## Bibliography auxiliary files (bibtex/biblatex/biber): +*.bbl +*.bcf +*.blg +*-blx.aux +*-blx.bib +*.run.xml + +## Build tool auxiliary files: +*.fdb_latexmk +*.synctex +*.synctex(busy) +*.synctex.gz +*.synctex.gz(busy) +*.pdfsync + +## Build tool directories for auxiliary files +# latexrun +latex.out/ + +## Auxiliary and intermediate files from other packages: +# algorithms +*.alg +*.loa + +# achemso +acs-*.bib + +# amsthm +*.thm + +# beamer +*.nav +*.pre +*.snm +*.vrb + +# changes +*.soc + +# comment +*.cut + +# cprotect +*.cpt + +# elsarticle (documentclass of Elsevier journals) +*.spl + +# endnotes +*.ent + +# fixme +*.lox + +# feynmf/feynmp +*.mf +*.mp +*.t[1-9] +*.t[1-9][0-9] +*.tfm + +#(r)(e)ledmac/(r)(e)ledpar +*.end +*.?end +*.[1-9] +*.[1-9][0-9] +*.[1-9][0-9][0-9] +*.[1-9]R +*.[1-9][0-9]R +*.[1-9][0-9][0-9]R +*.eledsec[1-9] +*.eledsec[1-9]R +*.eledsec[1-9][0-9] +*.eledsec[1-9][0-9]R +*.eledsec[1-9][0-9][0-9] +*.eledsec[1-9][0-9][0-9]R + +# glossaries +*.acn +*.acr +*.glg +*.glo +*.gls +*.glsdefs +*.lzo +*.lzs +*.slg +*.slo +*.sls + +# uncomment this for glossaries-extra (will ignore makeindex's style files!) +# *.ist + +# gnuplot +*.gnuplot +*.table + +# gnuplottex +*-gnuplottex-* + +# gregoriotex +*.gaux +*.glog +*.gtex + +# htlatex +*.4ct +*.4tc +*.idv +*.lg +*.trc +*.xref + +# hyperref +*.brf + +# knitr +*-concordance.tex +# TODO Uncomment the next line if you use knitr and want to ignore its generated tikz files +# *.tikz +*-tikzDictionary + +# listings +*.lol + +# luatexja-ruby +*.ltjruby + +# makeidx +*.idx +*.ilg +*.ind + +# minitoc +*.maf +*.mlf +*.mlt +*.mtc[0-9]* +*.slf[0-9]* +*.slt[0-9]* +*.stc[0-9]* + +# minted +_minted* +*.pyg + +# morewrites +*.mw + +# newpax +*.newpax + +# nomencl +*.nlg +*.nlo +*.nls + +# pax +*.pax + +# pdfpcnotes +*.pdfpc + +# sagetex +*.sagetex.sage +*.sagetex.py +*.sagetex.scmd + +# scrwfile +*.wrt + +# svg +svg-inkscape/ + +# sympy +*.sout +*.sympy +sympy-plots-for-*.tex/ + +# pdfcomment +*.upa +*.upb + +# pythontex +*.pytxcode +pythontex-files-*/ + +# tcolorbox +*.listing + +# thmtools +*.loe + +# TikZ & PGF +*.dpth +*.md5 +*.auxlock + +# titletoc +*.ptc + +# todonotes +*.tdo + +# vhistory +*.hst +*.ver + +# easy-todo +*.lod + +# xcolor +*.xcp + +# xmpincl +*.xmpi + +# xindy +*.xdy + +# xypic precompiled matrices and outlines +*.xyc +*.xyd + +# endfloat +*.ttt +*.fff + +# Latexian +TSWLatexianTemp* + +## Editors: +# WinEdt +*.bak +*.sav + +# Texpad +.texpadtmp + +# LyX +*.lyx~ + +# Kile +*.backup + +# gummi +.*.swp + +# KBibTeX +*~[0-9]* + +# TeXnicCenter +*.tps + +# auto folder when using emacs and auctex +./auto/* +*.el + +# expex forward references with \gathertags +*-tags.tex + +# standalone packages +*.sta + +# Makeindex log files +*.lpz + +# xwatermark package +*.xwm \ No newline at end of file diff --git a/doc/htway.html b/doc/htway.html new file mode 100644 index 0000000..6e72fda --- /dev/null +++ b/doc/htway.html @@ -0,0 +1,250 @@ + + + + + + Document + + + 
#!/usr/bin/env python
+         
+        # Python port of the HiTechnic HTWay for ev3dev
+        # Copyright (c) 2014-2015 G33kDude, David Lechner
+        # HiTechnic HTWay is Copyright (c) 2010 HiTechnic
+        import itertools
+        import os
+        import struct
+        import time
+        import pyudev
+         
+        # loop interval in seconds
+        WAIT_TIME = 0.008
+         
+        # ratio of wheel size compared to NXT 1.0
+        WHEEL_RATIO_NXT1 = 1.0 # 56mm
+        WHEEL_RATIO_NXT = 0.8 # 43.2mm (same as EV3)
+        WHEEL_RATIO_RCX = 1.4 # 81.6mm
+         
+        # These are the main four balance constants, only the gyro constants
+        # are relative to the wheel size. KPOS and KSPEED are self-relative to
+        # the wheel size.
+        KGYROANGLE = 7.5
+        KGYROSPEED = 1.15
+        KPOS = 0.07
+        KSPEED = 0.1
+         
+        # This constant aids in drive control. When the robot starts moving
+        # because of user control, this constant helps get the robot leaning in
+        # the right direction. Similarly, it helps bring robot to a stop when
+        # stopping.
+        KDRIVE = -0.02
+         
+        # Power differential used for steering based on difference of target
+        # steering and actual motor difference.
+        KSTEER = 0.25
+         
+        # If robot power is saturated (over +/- 100) for over this time limit
+        # then robot must have fallen.  In seconds.
+        TIME_FALL_LIMIT = 0.5
+         
+        # Gyro offset control
+        # The HiTechnic gyro sensor will drift with time.  This constant is
+        # used in a simple long term averaging to adjust for this drift.
+        # Every time through the loop, the current gyro sensor value is
+        # averaged into the gyro offset weighted according to this constant.
+        EMAOFFSET = 0.0005
+         
+        class Gyro:
+            @staticmethod
+            def get_gyro():
+                devices = list(pyudev.Context().list_devices(subsystem='lego-sensor') \
+                        .match_property("LEGO_DRIVER_NAME", EV3Gyro.DRIVER_NAME) \
+                        .match_property("LEGO_DRIVER_NAME", HTGyro.DRIVER_NAME))
+         
+                if not devices:
+                    raise Exception('Gyro not found')
+         
+                if devices[0].attributes.get('driver_name') == EV3Gyro.DRIVER_NAME:
+                    return EV3Gyro(devices[0])
+                elif devices[0].attributes.get('driver_name') == HTGyro.DRIVER_NAME:
+                    return HTGyro(devices[0])
+         
+            def __init__(self, device):
+                self.device = device
+                self.value0 = open(os.path.join(self.device.sys_path, 'value0'), 'r', 0)
+                self.offset = 0.0
+                self.angle = 0.0
+         
+            def calibrate(self):
+                self.angle = 0.0
+                total = 0
+                for i in range(10):
+                    total += self.get_rate()
+                    time.sleep(0.01)
+                average = total / 10.0
+                self.offset = average - 0.1
+         
+            def set_mode(self, value):
+                with open(os.path.join(self.device.sys_path, 'mode'), 'w') as f:
+                    f.write(str(value))
+         
+            def get_rate(self):
+                self.value0.seek(0)
+                return int(self.value0.read())
+         
+        class EV3Gyro(Gyro):
+            DRIVER_NAME = 'lego-ev3-gyro'
+         
+            def __init__(self, device):
+                Gyro.__init__(self, device)
+                self.set_mode("GYRO-RATE")
+         
+            def get_data(self, interval):
+                gyro_raw = self.get_rate()
+                self.offset = EMAOFFSET * gyro_raw + (1 - EMAOFFSET) * self.offset
+                speed = (gyro_raw - self.offset)
+                self.angle += speed * interval
+                return speed, self.angle
+         
+        class HTGyro(Gyro):
+            DRIVER_NAME = 'ht-nxt-gyro'
+         
+            def __init__(self, device):
+                Gyro.__init__(self, device)
+         
+            def get_data(self, interval):
+                gyro_raw = self.get_rate()
+                self.offset = EMAOFFSET * gyro_raw + (1 - EMAOFFSET) * self.offset
+                speed = (gyro_raw - self.offset) * 400.0 / 1953.0
+                self.angle += speed * interval
+                return speed, self.angle
+         
+        class EV3Motor:
+            def __init__(self, which=0):
+                devices = list(pyudev.Context().list_devices(subsystem='tacho-motor').match_attribute('driver_name', 'lego-ev3-l-motor'))
+         
+                if which >= len(devices):
+                    raise Exception("Motor not found")
+         
+                self.device = devices[which]
+                self.pos = open(os.path.join(self.device.sys_path, 'position'), 'r+',0)
+                self.duty_cycle_sp = open(os.path.join(self.device.sys_path,'duty_cycle_sp'), 'w', 0)
+                self.reset()
+         
+            def reset(self):
+                self.set_pos(0)
+                self.set_duty_cycle_sp(0)
+                self.send_command("run-direct")
+         
+            def get_pos(self):
+                self.pos.seek(0)
+                return int(self.pos.read())
+         
+            def set_pos(self, new_pos):
+                return self.pos.write(str(int(new_pos)))
+         
+            def set_duty_cycle_sp(self, new_duty_cycle_sp):
+                return self.duty_cycle_sp.write(str(int(new_duty_cycle_sp)))
+         
+            def send_command(self, new_mode):
+                with open(os.path.join(self.device.sys_path, 'command'),'w') as command:
+                    command.write(str(new_mode))
+         
+        class EV3Motors:
+            def __init__(self, left=0, right=1):
+                self.left = EV3Motor(left)
+                self.right = EV3Motor(right)
+                self.pos = 0.0
+                self.speed = 0.0
+                self.diff = 0
+                self.target_diff = 0
+                self.drive = 0
+                self.steer = 0
+                self.prev_sum = 0
+                self.prev_deltas = [0,0,0]
+         
+            def get_data(self, interval):
+                left_pos = self.left.get_pos()
+                right_pos = self.right.get_pos()
+         
+                pos_sum = right_pos + left_pos
+                self.diff = left_pos - right_pos
+         
+                delta = pos_sum - self.prev_sum
+                self.pos += delta
+         
+                self.speed = (delta + sum(self.prev_deltas)) / (4 * interval)
+         
+                self.prev_sum = pos_sum
+                self.prev_deltas = [delta] + self.prev_deltas[0:2]
+         
+                return self.speed, self.pos
+         
+            def steer_control(self, power, steering, interval):
+                self.target_diff += steering * interval
+                power_steer = KSTEER * (self.target_diff - self.diff)
+                power_left = max(-100, min(power + power_steer, 100))
+                power_right = max(-100, min(power - power_steer, 100))
+                return power_left, power_right
+         
+         
+        def main():
+            wheel_ratio = WHEEL_RATIO_NXT
+         
+            gyro = Gyro.get_gyro()
+            gyro.calibrate()
+         
+            print("balancing in ...")
+            for i in range(5, 0, -1):
+                print(i)
+                os.system("beep -f 440 -l 100")
+                time.sleep(1)
+            print(0)
+         
+            os.system("beep -f 440 -l 1000")
+            time.sleep(1)
+         
+            motors = EV3Motors()
+            start_time = time.time()
+            last_okay_time = start_time
+         
+            avg_interval = 0.0055
+         
+            for loop_count in itertools.count():
+                gyro_speed, gyro_angle = gyro.get_data(avg_interval)
+                motors_speed, motors_pos = motors.get_data(avg_interval)
+                print gyro_speed, gyro_angle, motors_speed, motors_pos
+                motors_pos -= motors.drive * avg_interval
+                motors.pos = motors_pos
+         
+                power = (KGYROSPEED * gyro_speed
+                        + KGYROANGLE * gyro_angle) \
+                        / wheel_ratio \
+                        + KPOS * motors_pos \
+                        + KDRIVE * motors.drive \
+                        + KSPEED * motors_speed
+         
+                left_power, right_power = motors.steer_control(power, 0, avg_interval)
+         
+                #print left_power, right_power
+         
+                motors.left.set_duty_cycle_sp(left_power)
+                motors.right.set_duty_cycle_sp(right_power)
+         
+                time.sleep(WAIT_TIME)
+         
+                now_time = time.time()
+                avg_interval = (now_time - start_time) / (loop_count+1)
+         
+                if abs(power) < 100:
+                    last_okay_time = now_time
+                elif now_time - last_okay_time > TIME_FALL_LIMIT:
+                    break
+         
+            motors.left.send_command('reset')
+            motors.right.send_command('reset')
+         
+        if __name__ == "__main__":
+            main()
+ + \ No newline at end of file diff --git a/doc/img/ang_run_1.png b/doc/img/ang_run_1.png new file mode 100644 index 0000000..ea4c353 Binary files /dev/null and b/doc/img/ang_run_1.png differ diff --git a/doc/img/ang_sim_1.png b/doc/img/ang_sim_1.png new file mode 100644 index 0000000..a7a144c Binary files /dev/null and b/doc/img/ang_sim_1.png differ diff --git a/doc/img/cart_2.png b/doc/img/cart_2.png new file mode 100644 index 0000000..4344eb9 Binary files /dev/null and b/doc/img/cart_2.png differ diff --git a/doc/img/pendulum_2.png b/doc/img/pendulum_2.png new file mode 100644 index 0000000..eaafd86 Binary files /dev/null and b/doc/img/pendulum_2.png differ diff --git a/doc/img/run1.png b/doc/img/run1.png new file mode 100644 index 0000000..43f4ef2 Binary files /dev/null and b/doc/img/run1.png differ diff --git a/doc/img/run2.png b/doc/img/run2.png new file mode 100644 index 0000000..93aa2b4 Binary files /dev/null and b/doc/img/run2.png differ diff --git a/doc/img/run3.png b/doc/img/run3.png new file mode 100644 index 0000000..7855202 Binary files /dev/null and b/doc/img/run3.png differ diff --git a/doc/img/sim_run_alllimits.png b/doc/img/sim_run_alllimits.png new file mode 100644 index 0000000..82c6d6a Binary files /dev/null and b/doc/img/sim_run_alllimits.png differ diff --git a/doc/img/sim_run_alllimits2.png b/doc/img/sim_run_alllimits2.png new file mode 100644 index 0000000..5240fe3 Binary files /dev/null and b/doc/img/sim_run_alllimits2.png differ diff --git a/doc/img/sim_run_pidcontroldelay.png b/doc/img/sim_run_pidcontroldelay.png new file mode 100644 index 0000000..ebfa6d8 Binary files /dev/null and b/doc/img/sim_run_pidcontroldelay.png differ diff --git a/doc/img/sim_run_resolutionlimit1.png b/doc/img/sim_run_resolutionlimit1.png new file mode 100644 index 0000000..9455c25 Binary files /dev/null and b/doc/img/sim_run_resolutionlimit1.png differ diff --git a/doc/img/sim_run_speedlimit.png b/doc/img/sim_run_speedlimit.png new file mode 100644 index 0000000..58bb30e Binary files /dev/null and b/doc/img/sim_run_speedlimit.png differ diff --git a/doc/img/sim_run_speedlimit_ms2.png b/doc/img/sim_run_speedlimit_ms2.png new file mode 100644 index 0000000..9ad72d4 Binary files /dev/null and b/doc/img/sim_run_speedlimit_ms2.png differ diff --git a/doc/pendulum.bib b/doc/pendulum.bib new file mode 100644 index 0000000..fef138c --- /dev/null +++ b/doc/pendulum.bib @@ -0,0 +1,32 @@ +@MISC{instruction, + author = "Arne van Iterson", + title = "Building instructions", + url = {https://arnweb.nl/gitea/arne/EV5_Modcon/releases/download/pre-alpha/Instruction.pdf} +}, + +@misc{htway, + author = {David Lechner}, + howpublished = {Github}, + title = {Python port of the HiTechnic HTWay for ev3dev}, + year = {2010}, + url = {https://gist.github.com/dlech/11098915}, +} + +@article{motorcomp, + author = {Philippe Hurbain}, + title = {LEGO® 9V Technic Motors compared characteristics}, + year = {2012}, + url = {https://www.philohome.com/motors/motorcomp.htm} +} + +@ARTICLE{writers, + author = "Los Angles Times", + title = "Writers’ strike: What happened, how it ended and its impact on Hollywood", + url = {https://www.latimes.com/entertainment-arts/business/story/2023-05-01/writers-strike-what-to-know-wga-guild-hollywood-productions} +}, + +@ARTICLE{support, + author = "Washington Post", + title = "ChatGPT provided better customer service than his staff. He fired them.", + url = {https://www.washingtonpost.com/technology/2023/10/03/ai-customer-service-jobs/} +} \ No newline at end of file diff --git a/doc/pendulum.tex b/doc/pendulum.tex index 255ac1a..5f66f89 100644 --- a/doc/pendulum.tex +++ b/doc/pendulum.tex @@ -2,16 +2,20 @@ \documentclass{article} % \documentclass{} is the first command in any LaTeX code. It is used to define what kind of document you are creating such as an article or a book, and begins the document preamble +\usepackage[english]{babel} \usepackage[a4paper,top=2cm,bottom=2cm,left=2cm,right=2cm,marginparwidth=1.75cm]{geometry} -\usepackage{amsmath} % \usepackage is a command that allows you to add functionality to your LaTeX code +\usepackage{amsmath} \usepackage{multicol} \usepackage{textcomp} \usepackage{graphicx} \usepackage{float} +\usepackage{lipsum} +\usepackage{hyperref} +\usepackage{listings} \title{ - Single inverted pendulum\\~\ - \large{Modeling and Control course} + Single inverted pendulum\linebreak + \large{EV5 - Modelling and Control} } \author{Arne van Iterson} % Sets authors name @@ -23,47 +27,886 @@ \begin{abstract} - Very interesting Introduction - + Electrical engineering students at the University of applied sciences Utrecht are assigned to simulate, build and control an inherently unstable system during the 'Modelling and Control' course of the third year. This paper describes the process of doing just that for a single inverted pendulum. The pendulum itself is built using LEGO\textsuperscript{\tiny\textregistered} Mindstorms EV3. The control system is a PID-loop is implemented using Python. A digital twin, also built in Python, is used to simulate the system. \end{abstract} +\noindent\makebox[\linewidth]{\rule{\linewidth}{0.4pt}} + \begin{multicols}{2} % creates two columns % Theory \section{Theory} - The system that will be discussed in this paper is a single inverted pendulum. Common implementations of this system include a cart, this is not the case in this paper. + \label{section:theory} + The system that will be discussed in this paper is a single inverted pendulum. Common implementations of this system include a cart of some sort that carries the pendulum on a pivot point over a defined path. + + While the system could theoretically stand up straight and never fall over, in practice the system is inherently unstable and will fall over if not controlled. The system can be controlled using a motor that moves the cart up and down the path to compensate for the pendulums falling movement. This method of control only requires a sensor to measure the angle of the pendulum. - \begin{figure}[H] - \includegraphics[width=\linewidth]{../res/segway.jpg} - \caption{Segway\textsuperscript{\tiny\textregistered} in use} - \label{fig:segway} - \end{figure} + \begin{figure}[H] + \includegraphics[width=\linewidth]{../res/drawio/system_diagram.png} + \caption{System with (left) and without cart (right)} + \label{fig:sysdiag} + \end{figure} + + The system used in this case of this assignment is slightly different; Instead of a cart that carries the pendulum, the pendulum itself is equipped with a motor and the angle is measured using a gyroscope. + This difference has been visualised in figure \ref{fig:sysdiag}. The physics of the system are practically identical, however, it allows for some minor simplifications which will be discussed in section \ref{section:model}. + + A practical example of this system is a Segway\textsuperscript{\tiny\textregistered} People Transporter (figure \ref{fig:segway}). + + \begin{figure}[H] + \includegraphics[width=\linewidth]{../res/segway.jpg} + \caption{Two Segway\textsuperscript{\tiny\textregistered} PT units in use} + \label{fig:segway} + \end{figure} - A practical example of this system is a Segway\textsuperscript{\tiny\textregistered} (figure \ref{fig:segway}). The system is inherently unstable, meaning that if the pendulum is not controlled, it will fall over. + \section{Model} + \label{section:model} + The physics behind the system in question have been described using the Lagrangian method. The Lagrangian method is a way of describing the dynamics of a system using the kinetic and potential energy of the system. + + \subsection{Lagrangian} + Usually the model for this system is divided into two parts: The cart and the pendulum. The cart is usually considered to be a mass $M$ that moves horizontally along a track. The pendulum is usually considered to be a mass $m$ that is attached to the cart using a pivot point. In this case however, since the motor is attached to the pendulum and the pivot point is the driving shaft of the motor itself, there is no separate mass for the cart. This means that $M$ is equal to $m$ and the formula can be simplified slightly. + + \newpage + + The kinetic energy of the system is made up from by the $x$ movement of the system and the falling of the pendulum. Since the system can only move in the $x$ direction using the motor, the kinetic energy is simply: + \begin{equation} + T_{1} = \frac{1}{2} m v^2 = \frac{1}{2} m \dot x^2 + \end{equation} + Where $x$ is the displacement of the system in the $x$ direction and $m$ is the mass of the system. \\ + + The kinetic energy of the pendulum is influenced by both the $x$ and $y$ movements of the pendulum and therefore: + \begin{equation*} + T_{2} = \frac{1}{2} m (\dot x_p^2 + \dot y_p^2) + \end{equation*} + The $x$ and $y$ position of the top of the pendulum can be calculated using the length and angle, this assumes that the pendulum is upright at $\frac{1}{2}\pi$: + \begin{equation*} + \begin{split} + x_{p} & = x + l\cdot sin(\theta) \\ + y_{p} & = -l\cdot cos(\theta) + \end{split} + \end{equation*} + The $x$ and $y$ velocity can be derived from the position: + \begin{equation*} + \begin{split} + \dot{x}_{p} & = \dot{x} + l\cdot cos(\theta)\cdot \dot{\theta} \\ + \dot{y}_{p} & = l\cdot sin(\theta)\cdot \dot{\theta} + \end{split} + \end{equation*} + Putting this all together gives the kinetic energy of the pendulum: + \begin{align*} + T_{2} &= \frac{1}{2} m (\dot x_s^2 + \dot y_s^2) \\ + T_{2} &= \frac{1}{2} m((\dot x + \ell \dot \theta cos(\theta))^2 + (\ell \dot \theta sin(\theta))^2) \\ + T_{2} &= \frac{1}{2} m(\dot x^2 + 2 \dot x \ell \dot \theta cos(\theta) + \ell^2 \dot \theta^2 cos(\theta)^2 + \ell^2 \dot \theta^2 sin(\theta)^2) \\ + T_{2} &= \frac{1}{2} m(\dot x^2 + 2 \dot x \ell \dot \theta cos(\theta) + \ell^2 \dot \theta^2(cos(\theta)^2 + sin(\theta)^2)) \\ + \intertext{Given that $cos(x)^2 + sin(x)^2 = 1$:} + T_{2} &= \frac{1}{2} m(\dot x^2 + 2 \dot x \ell \dot \theta cos(\theta) + \ell^2 \dot \theta^2) + \intertext{Combining $T_1$ and $T_2$ gives the total kinetic energy of the system:} + T & = T_1 + T_2 \\ + T & = \frac{1}{2} m \dot x^2 + \frac{1}{2} m(\dot x^2 + 2 \dot x \ell \dot \theta cos(\theta) + \ell^2 \dot \theta^2) \\ + T & = m \dot x^2 + \frac{1}{2} m(2 \dot x \ell \dot \theta cos(\theta) + \ell^2 \dot \theta^2) \\ + T & = m \dot x^2 + m \dot x \ell \dot \theta cos(\theta) + \frac{1}{2} m\ell^2 \dot \theta^2 + \end{align*} + The potential energy of the system is: + \begin{equation} + V = -mg\ell cos(\theta) + \end{equation} + This makes the full Lagrangian of the system: + \begin{equation} + \begin{split} + \mathcal{L} & = T - V \\ + \mathcal{L} & = m\dot x^2 + m \dot x \ell \dot \theta cos(\theta) + \frac{1}{2}m \ell^2 \dot \theta^2 + mg\ell cos(\theta) + \end{split} + \end{equation} + + Since the falling pendulum influences the moving pivot point ("cart") and vice versa, the Lagrangian will have to be solved for both the angular acceleration $\ddot\theta$ of the pendulum and the horizontal acceleration of the pivot $\ddot{x}$ + + \subsection{Horizontal acceleration} + The horizontal acceleration of the pivot point can be derived from the Lagrangian using the following formula: + \begin{equation} + \label{eq:horac} + \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot x} \right) = \frac{\partial \mathcal{L}}{\partial x} + \end{equation} + Filling in for $\dot x$ + % Dependant on xdot + \begin{equation} + \label{eq:horacxdot} + \begin{split} + &\ m\dot x^2 + m \dot x \ell \dot \theta cos(\theta) \\ + \frac{\partial \mathcal{L}}{\partial \dot x} =&\ 2m\dot x + m \ell \dot \theta cos(\theta) \\ + \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot x} \right) =&\ 2m\ddot x + m \ell \ddot \theta cos(\theta) - m \ell \dot \theta sin(\theta) + \end{split} + \end{equation} + Filling in for $x$ + % Dependant on x + \begin{equation} + \label{eq:horacx} + \frac{\partial \mathcal{L}}{\partial x} = 0 + \end{equation} + Combining formulas \ref{eq:horacxdot} and \ref{eq:horacx} using formula \ref{eq:horac} gives: + \begin{equation} + 2m\ddot x + m \ell \ddot \theta cos(\theta) - m \ell \dot \theta sin(\theta) = 0 + \end{equation} + Which gives: + \begin{equation} + \begin{split} + 2m\ddot x &= - m \ell \ddot \theta cos(\theta) + m \ell \dot \theta sin(\theta) \\ + \ddot x &= -\frac{m \ell \ddot \theta cos(\theta) + m \ell \dot \theta sin(\theta)}{2m} + \end{split} + \end{equation} + + \subsection{Angular acceleration} + \begin{equation} + \label{eq:angac} + \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot \theta} \right) = \frac{\partial \mathcal{L}}{\partial \theta} + \end{equation} + Filling in for $\dot \theta$ + % Dependant on thetadot + \begin{equation} + \label{eq:angacthetadot} + \begin{split} + &\ m \dot x \ell \dot \theta cos(\theta) + \frac{1}{2}m \ell^2 \dot \theta^2 \\ + \frac{\partial \mathcal{L}}{\partial \dot \theta} =&\ m \dot x \ell cos(\theta) + m \ell^2 \dot\theta\\ + \frac{d}{dt} \left(\frac{\partial \mathcal{L}}{\partial \dot \theta} \right) =&\ m \ddot x \ell cos(\theta) - m \dot x \ell sin(\theta) + m \ell^2 \ddot\theta\\ + \end{split} + \end{equation} + Filling in for $\theta$ + % Dependant on theta + \begin{equation} + \label{eq:angactheta} + \begin{split} + &\ m \dot x \ell \dot \theta cos(\theta) + mg\ell cos(\theta) \\ + \frac{\partial \mathcal{L}}{\partial \theta} =&\ -m \dot x \ell \dot \theta sin(\theta) - mg\ell sin(\theta) \\ + \end{split} + \end{equation} + Combining formulas \ref{eq:angacthetadot} and \ref{eq:angactheta} using formula \ref{eq:angac} gives: + \begin{equation} + \begin{split} + m \ell^2 \ddot\theta =& - m \ddot x \ell cos(\theta) + m \dot x \ell sin(\theta) -m \dot x \ell \dot \theta sin(\theta) - mg\ell sin(\theta) \\ + \ell \ddot\theta =& - \ddot x cos(\theta) + \dot x sin(\theta) -\dot x \dot \theta sin(\theta) - g sin(\theta) \\ + \ddot\theta =& - \frac{\ddot x cos(\theta) - \dot x sin(\theta) + \dot x \dot \theta sin(\theta) + g sin}{l} + \end{split} + \end{equation} + + % Make clear that the simulation is built from the ground up and pygame is only used for the UI + \subsection{Simulation} + The formula's in the previous sections have then been put into a purpose built Python simulation suite called Pendulum Simulator 4000. The source code for this can be found in the appendix and on git. + \begin{figure}[H] + \includegraphics[width=\linewidth]{../res/screenshot2.png} + \caption{Main window of Pendulum Simulator 4000} + \label{fig:simmain} + \end{figure} + % Hardware setup \section{Setup} - To simplify hardware design, the system will be built and run on a LEGO\textsuperscript{\tiny\textregistered} Mindstorms EV3 development kit. The programmable brick itself, which will be referred to as "EV3", runs a Linux distribution known as ev3dev. Ev3dev allows for various programming languages to be used including Python and C. Both of which will be further explained later. + To simplify hardware design, the system will be built and run on a LEGO\textsuperscript{\tiny\textregistered} Mindstorms EV3 development kit. The programmable brick itself, hereinafter referred to as "EV3", runs a Linux distribution known as ev3dev. Ev3dev allows for various programming languages to be used including Python and C. Both of which will be further explained in section \ref{section:control}. - To reproduce the robot in Figure \ref{fig:balancer}, please refer to the building instructions in the appendix. + To reproduce the robot in Figure \ref{fig:balancer}, please refer to the building instructions. \cite{instruction}. - \begin{figure}[H] - \includegraphics[width=\linewidth]{../res/bricklink/modcon_balancer2.png} - \caption{Balancing robot} - \label{fig:balancer} - \end{figure} + \begin{figure}[H] + \includegraphics[width=\linewidth]{../res/bricklink/modcon_balancer2.png} + \caption{Balancing robot} + \label{fig:balancer} + \end{figure} -\section{Model} -\subsection{Linearisation} + \section{Control} + \label{section:control} + In the simulation it was found that the system was best controlled using a PD-loop, this was to be expected since the system includes integration by it's very nature. The control loop was implemented using Python and the ev3dev library. The source code for this can be found in the appendix. -\section{Control} -dsaf + The choice of hardware has significantly influenced the control system. Unfortunately, it has brought along a lot of issues with the physical model. The gyroscope used has a very limited resolution of 1 degree (1 deg) for angle and 1 degree per second (1 deg/sec) for angular velocity. Add the fact that the sensor drifts about 0,5 deg/sec and the system becomes very unstable very quickly. -\section{Validation} -fsafdfa + This issue would probably have been overcome if the system was easier to program, in its current state the EV3 runs a limited version of Python that takes about two minutes to load every time the program is run. Adding any reasonable logging, through the file system or through dumping python arrays, makes the timing too slow to keep the control loop running fast enough causing the robot to lose balance fairly quickly. + + There has been some success in using the C programming language and a compiling toolchain has been set-up, however, the control of simple features like the motor is very low level. The plan was to first get it working in Python and then port it to C but time has unfortunately run out. -\section{Conclusion} -asfjdklasdfas + In order to get the system to balance at all, the control loop had to be altered to include the position and speed of the motors as has been done in the HTWay project by D. Lechner \cite{htway}. The theory behind this is that if the system is still falling in a certain direction while the motors are already moving to compensate that, the movement speed should be increased. In a way the system now uses a cascaded control loop to balance. + + A small set of data points has been collected by pre-assigning memory to a numpy array and writing to it directly, then dumping this to non-volatile memory. This gives some insight into the issues with the system, however, this also limits the run time of the system. + + \section{Validation} + The behaviour of the system is as follows: + + \begin{figure}[H] + \includegraphics[width=\linewidth]{./img/ang_run_1.png} + \caption{Run 1} + \end{figure} + + The fact that the system oscillates is a typical sign of an ill-tuned control loop. The PID values have been tuned to the best of my ability but the system is still very unstable. Tuning the system has proven difficult due to the inclusion of the motor speed and position in the control loop. + + The system is able to say upright using the HTWay project by D. Lechner \cite{htway}. Some of the control logic has been used in this project, however the HTWay project includes some filtering logic for the gyro sensor that I did not yet understood well enough to implement. + \\\\ + The fact that I was unable to get the system to balance properly should not change the fact that the simulation should behave the same given the same characteristics and control constants as the physical model, which is not the case. + + Running the simulation tuned exactly the same as the physical model reveals an interesting issue; The simulated system is impossible to topple over, which is clearly not the case with the physical model. + \\\\ + The simulation and the model differ in a few ways: + \begin{itemize} + \item The simulation does not include sensor drift or resolution + \item The simulation does not include the physical limitations of the motor + \item The simulation loop runs significantly faster than the physical model does + \end{itemize} + + Drawing any comparisons between the two will be largely useless at this point, so instead we will limit the simulation to check if the simulation will then mimic the actual behaviour. + + \subsection{Limiting the simulated system} + First, the resolution of the gyroscope will be limited to 1 degree. Converting this to radians gives a resolution of 0.0175 rad. The result can be seen in \ref{fig:simreslim1}. Note that the graph has been cropped, showing only a stable part of the simulation. + \begin{figure}[H] + \includegraphics[width=\linewidth]{./img/sim_run_resolutionlimit1.png} + \caption{Simulation run with angle resolution limit} + \label{fig:simreslim1} + \end{figure} + + The limited resolution does cause the system to oscillate slightly, but not escalating as expected. Additionally, the system is still impossible to topple over. + \\\\ + Next, the simulation will be limited in update speed. This means that the physics simulation will continue while the control loop can only update every certain interval. The log file from the physical model reveals that the average update time for the control loop is about 21 ms or 47 frames per second (fps). + + Currently, the simulation runs at 60 fps, the update call for the pendulum control is called every frame. The pendulum control has been changed to only actually control the system after the dt value is more or equal to 21 ms. Because the program runs at 60 fps or roughly 16 ms per frame, the actual PID loop will be called every 32 ms; This is longer than required, but should demonstrate the effect of the delay. The result can be seen in \ref{fig:simpidlim1}. + \begin{figure}[H] + \includegraphics[width=\linewidth]{./img/sim_run_pidcontroldelay.png} + \caption{Simulation run with PID control delay} + \label{fig:simpidlim1} + \end{figure} + + Interestingly, the limited update frequency of the PID controller seems to act as a low-pass filter, reducing the oscillations in the system. The user is now able to topple the system over. + \\\\ + Next, the simulation will be limited in the speed at which it can change the horizontal position of the pendulum. The motor speed on the EV3 is measured in an arbitrary value between -1050 and 1050. The max speed of the motor is about 175 rpm as stated in the research by P. Hurbain. \cite{motorcomp}. Combined with the 44 mm diameter of the used wheels gives the maximum speed of the robot as: + \begin{equation} + \begin{split} + r_{wheel} &= 22 \text{[mm]}\\ + s_{circumference} &= 2\pi r_{wheel} = 138.23 \text{mm} = 0.1382 \text{m}\\ + v_{max} &= (175 / 60) \cdot 0.1382 \\ + &= 2.92 \cdot 0.1382 = 0.403 \text{m/s} + \end{split} + \end{equation} + + The simulation has been changed to never exceed the speed limit of 0.403 m/s even if the controller asks for it, just like the EV3 motors would do. Limiting the simulation to 0.403 m/s did little to change the behaviour of the simulation. As can be seen in \ref{fig:simslms}, the system never reaches this speed in the first place. + + \begin{figure}[H] + \includegraphics[width=\linewidth]{./img/sim_run_speedlimit.png} + \caption{Simulation run with speed limit} + \label{fig:simslms} + \end{figure} + + Interesting however is the acceleration as can be seen in \ref{fig:simslms2}. + + \begin{figure}[H] + \includegraphics[width=\linewidth]{./img/sim_run_speedlimit_ms2.png} + \caption{Simulation run with speed limit} + \label{fig:simslms2} + \end{figure} + + The initial kick in acceleration is the attempt to topple the system over by the user. The response is quite harsh, reaching up to 7 m/s\textsuperscript{2} which is only slightly impossible for the physical system to achieve. Unfortunately, this is also quite difficult to limit, since the acceleration of the motor would require a measuring setup to determine. + + If we make a generous assumption that the motor can reach it's top speed within 0.15 seconds, the acceleration would be: + \begin{equation} + \begin{split} + v_{max} &= 0.403 \text{m/s} \\ + t_{max} &= 0.1 \text{s} \\ + a_{max} &= \frac{v_{max}}{t_{max}} = \frac{0.403}{0.15} = 2.69\text{m/s}^2 + \end{split} + \end{equation} + + This does cause the system to fall over in the simulation as expected, though not as quickly as the physical model and with less oscillations. + + \begin{figure}[H] + \includegraphics[width=\linewidth]{./img/sim_run_alllimits.png} + \caption{Simulation run with all limits in place} + \label{fig:simalllim1} + \end{figure} + + Lastly, the sensor will be offset with about 0,5 degrees per second which will be updated every time the control loop is run. The results can be seen in \ref{fig:simdrift1}. + + \begin{figure}[H] + \includegraphics[width=\linewidth]{./img/sim_run_alllimits2.png} + \caption{Simulation run with sensor drift} + \label{fig:simdrift1} + \end{figure} + + With that, we have successfully created a digital twin for a poorly tuned system. The system is now no longer able to balance properly, it oscillates and falls over just like the physical model. + + \section{Conclusion} + Due to issues with the physical model, this entire project has run slightly backwards. Having to alter the control loop significantly nearing the end of the project has caused the physical model to behave differently than the simulation to a point where comparing them was not really possible. + The simulation has been altered to include the limitations of the physical model as closely as possible, which has been successful and the resulting behaviour is comparable to the physical model. So while the system does not balance, at least it makes sense. + + \section{Recommendations} + \begin{enumerate} + \item Do not make a robot using LEGO. + \end{enumerate} \end {multicols} +\newpage + +\section{Ethics of Artificial Intelligence} +\label{appendix:ethics} + In Utilitarian ethics the greatest good for the greatest number is the most important. Whatever happens along the way does not matter as long as the outcome is good for the most people or society as a whole. The last couple of years have seen great development in the field of Artificial Intelligence, which brings along the question of how and if this technology should be used from a utilitarian perspective. + \\\\ + The answer to this question is not as simple as it may seem. On the one hand, AI can be used to automate tasks that are dangerous or tedious for humans. This would allow them to focus on more important tasks and would reduce the amount of work related injuries; Which would be favourable according to utilitarian ethics. This also means however, that it would replace work that is currently done by humans, therefore leading to a loss of jobs and, potentially, a group of people that are unable to find work. + \\\\ + Then again, we are currently experiencing the development of AI, which means that what is happening now might just be part of the AI development that will eventually lead to the greatest thing humans have ever done that will eventually benefit everyone. In that case, we might just have to accept the loss of jobs and the other negative effects of AI development as a necessary evil until the great-for-everyone endpoint is reached. + \\\\ + This great-for-everyone endpoint however, is not a certainty and if it is ever reached, it might take a very long time. Therefore, I am of the opinion that not all types of AI should be treated with the same ethical standards. + \\\\ + In the context of modelling and control, I think most of the work described in this paper, maybe with the exception of building the physical model itself, will no longer be done by a person in the near future. Writing this paper and parts of the software in VSCode with Github Copilot Enabled has taught me that with a more limited context compared to generic AI, like ChatGPT, the AI 'understood' more complex code structures and physics than I expected. At the moment, I don't think it would be a great idea to run Copilot code without any checking by a human, but I do think that it will be possible in the near future. + \\\\ + There is a lot more to engineering than just coding, therefore I don't think engineers as a whole will be replaced at all. AI development towards the aid of developers will likely not directly cost any jobs. There have been cases where it has or likely will, especially creative jobs, journalism and administrative jobs. \cite{writers} This type of development, I believe, should not stop but should be kept in check. We have seen the writers strikes and the companies firing their entire support staff after OpenAI opened the GPT-3 API. \cite{support} While it might be part of the upcoming great-for-everyone endpoint, right now I do not believe we should allow AI to completely replace jobs with half-baked solutions for the sake of profit. + \\\\ + +\newpage + +\bibliographystyle{IEEEtran} +\bibliography{pendulum} + +\newpage + +\appendix + \section{Simulation source code} + \subsection{sim.py} + \begin{lstlisting}[language=Python] + # Pendulum simulator 4000 + # Arne van Iterson, 2023 + + # Imports + import pygame_widgets + import pygame + from pygame_widgets.slider import Slider + from pygame.math import Vector2 + import math + + # pygame setup + pygame.init() + screen = pygame.display.set_mode((1280, 720)) + clock = pygame.time.Clock() + running = True + update = True + pole = Vector2(screen.get_rect().center) # center of screen + + # Own objects must be imported after pygame init + from pendulum import Pendulum + from uiHelpers import * + + # UI helpers + ui = SimUI(screen, pole) + + # Pendulum setup + # Start angle in radians, length, mass, color + pendulum = Pendulum(0, 2, 0, 0.25, "red") + pendulum.reset() + dx = 0 # x offset + dt = 1 # delta time + + # Sliders + slider_kp = Slider(screen, 910, 590, 320, 16, initial=pendulum.kp, min=0, max=2, step=0.001) + slider_ki = Slider(screen, 910, 620, 320, 16, initial=pendulum.ki, min=0, max=0.5, step=0.001) + slider_kd = Slider(screen, 910, 650, 320, 16, initial=pendulum.kd, min=0, max=0.5, step=0.001) + + # Gametime + rt = 10 # run time + highscore = 0 + + + # Metadata values + def meta(): + ui.meta(pendulum.theta[pendulum.index], "Theta") + ui.meta(pendulum.a_ang[pendulum.index], "Angular acceleration") + ui.meta(pendulum.dx, "dx") + ui.meta(pendulum.a_cart[pendulum.index], "Cart acceleration") + + ui.meta(pendulum.pid, "Control") + + ui.meta(pendulum.kd, "Kd") + ui.meta(pendulum.ki, "Ki") + ui.meta(pendulum.kp, "Kp") + + ui.meta(dt, "dt") + + ui.meta(not update, "Paused") + ui.meta(rt / 1000, "Run time [s]") + ui.meta(highscore / 1000, "Highscore [s]") + + + while running: + events = pygame.event.get() + ### User controls ### + for event in events: + # Quit + if event.type == pygame.QUIT: + running = False + elif event.type == pygame.KEYDOWN: + # Quit + if event.key == pygame.K_ESCAPE: + running = False + # Reset simulation + elif event.key == pygame.K_SPACE: + pendulum.reset() + rt = 0 + # Pause simulation + elif event.key == pygame.K_p: + update = not update + # Display plot if simulation is not running + elif event.key == pygame.K_g: + if pendulum.fallen: + pendulum.plot() + else: + update = False + pendulum.plot() + # Toggle PID controller + elif event.key == pygame.K_c: + pendulum.pid = not pendulum.pid + + # Move pendulum + keys = pygame.key.get_pressed() + if keys[pygame.K_LEFT] or keys[pygame.K_a]: + pendulum.a_cart[pendulum.index] -= 4 + if keys[pygame.K_RIGHT] or keys[pygame.K_d]: + pendulum.a_cart[pendulum.index] += 4 + + # Draw grid + ui.grid(50, 0, 15) + + # Update PID values + pendulum.kp = slider_kp.getValue() + pendulum.ki = slider_ki.getValue() + pendulum.kd = slider_kd.getValue() + + # Update pendulum + if not pendulum.fallen: + if update: + rt += dt + pendulum.update(dt) + else: + ui.gameover(rt) + + # Update highscore + if rt > highscore: + highscore = rt + + # Draw metadata + ui.update(dt) + meta() + + # Draw pendulum + dx = (pendulum.dx, 0) + pygame.draw.line(screen, pendulum.color, pole + dx, pole + pendulum.vector + dx, 3) + pygame.draw.circle(screen, "black", pole + dx, 15, 3) + + # Draw frame + pygame_widgets.update(events) + pygame.display.flip() + dt = clock.tick(60) # limits FPS to 120 + + pygame.quit() + + \end{lstlisting} + \subsection{pendulum.py} + \begin{lstlisting}[language=Python] + from pygame.math import Vector2 + import math + import numpy as np + import random + import matplotlib.pyplot as plt + + # Constants + C_GRAVITY = 9.81 # m/s^2 + C_MTPRATIO = 100 # Pixels per meter + C_P_ANG_START = 1 / 1000 * math.pi + C_FALL_ANG = 52.5 / 100 * math.pi + + + class Pendulum: + def __init__(self, theta, length, dx, mass, color): + """ + Initialize a Pendulum object. + + Parameters: + theta (float): Angle [rad]. + length (float): Length of the pendulum [m]. + dx (float): Horizontal displacement of the "cart" from the center [m]. + mass (float): Mass of the pendulum for physics calculations [kg]. + color (str): Display color. + + Returns: + None + """ + ## Game variables + self.vector = None # Vector2 object + self.fallen = False # Stop when pendulum falls over + + ## Physics variables + self.index = 0 # Index helper for plotting graphs + self.theta = [theta] # Angle in radians + self.a_ang = [0] # Angular acceleration + self.v_ang = [0] # Angular velocity + + self.dx = dx # Horizontal displacement of "cart" from center + self.a_cart = [0] # Acceleration of cart + self.v_cart = [0] # Velocity of cart + self.s_cart = [0] # Displacement of cart [m] + + # self.r_factor = 0.50 # Damping factor + + self.length = length # Length of pendulum + self.mass = mass # Mass of pendulum for physics + self.color = color # Display color + + ## PID variables + self.pid = True + self.kp = 1.3 + self.ki = 0.0 + self.kd = 0.1 + + def update(self, dt): + self.doMath(dt) + self.vector = Vector2.from_polar( + ((self.length * C_MTPRATIO), math.degrees(self.theta[self.index] - (0.5 * math.pi))) + ) + + if self.pid: + self.pidControl(dt) + + if abs(self.theta[self.index]) == C_FALL_ANG: + self.fallen = True + + # def update(self, dt): + # """ + # Update the pendulum's state based on the elapsed time. + + # Parameters: + # - dt (float): The elapsed time in milliseconds. + + # Returns: + # None + # """ + # a_ang = (-(C_GRAVITY * math.sin(self.theta)) / (self.length)) - (self.r_factor * self.v_ang) # Angular acceleration + + # v_ang = a_ang * (dt/1000) + self.v_ang # Integrate acceleration to get velocity + # s_ang = v_ang * (dt/1000) # Angular displacement + + # self.theta += s_ang # Update value + + # self.vector = Vector2.from_polar(((self.length * 150), math.degrees(self.theta + math.pi/2))) + + # self.a_ang = a_ang # Update value + # self.v_ang = v_ang # Update value + + def doMath(self, dt): + ### ANGLE ### + ang_term1 = self.a_cart[self.index] * math.cos(self.theta[self.index]) + ang_term2 = self.v_cart[self.index] * math.sin(self.theta[self.index]) + ang_term3 = ( + self.v_cart[self.index] # Previous cart velocity + * self.v_ang[self.index] # previous angle velocity + * math.sin(self.theta[self.index]) # Sin previous angle + ) + ang_term4 = C_GRAVITY * math.sin(self.theta[self.index]) + + # Angular acceleration + self.a_ang.append( + (ang_term1 - ang_term2 + ang_term3 - ang_term4) / -(self.length) + ) + + # Integrate acceleration to get velocity + self.v_ang.append( + self.v_ang[self.index] # Previous velocity + + (self.a_ang[self.index + 1] * (dt / 1000)) + ) + + # Angular displacement + self.theta.append( + self.theta[self.index] # Previous angle + + (self.v_ang[self.index + 1] * (dt / 1000)) + ) + + # Limit fall of pendulum + self.theta[self.index + 1] = self.clamp( + self.theta[self.index + 1], -C_FALL_ANG, C_FALL_ANG + ) + + ### CART ### + cart_term1 = ( + self.mass # Mass + * self.length # Length + * self.a_ang[self.index + 1] # Current angle acceleration + * math.cos(self.theta[self.index + 1]) # Current angle + ) + cart_term2 = ( + self.mass # Mass + * self.length # Length + * self.v_ang[self.index + 1] # Current angle velocity + * math.sin(self.theta[self.index + 1]) # Current angle + ) + + # Cart acceleration + self.a_cart.append((-cart_term1 + cart_term2) / (2 * self.mass)) + + # Integrate acceleration to get velocity + self.v_cart.append( + self.v_cart[self.index] # Previous velocity + + (self.a_cart[self.index + 1] * (dt / 1000)) + ) + + # Cart displacement + self.s_cart.append( + self.s_cart[self.index] # Previous displacement + + (self.v_cart[self.index + 1] * (dt / 1000)) + ) + self.dx = self.s_cart[self.index + 1] * C_MTPRATIO # Convert to pixels + + # Update index + self.index += 1 + + def clamp(self, n, minn, maxn): + return max(min(maxn, n), minn) + + def reset(self): + self.index = 0 + + self.a_ang = [0] + self.v_ang = [0] + self.dx = [0] + self.theta = [random.choice([1, -1]) * C_P_ANG_START] + + self.a_cart = [0] + self.v_cart = [0] + self.s_cart = [0] + self.fallen = False + self.update(0) + + def plot(self): + fig, axs = plt.subplots(2, 2) + fig.suptitle("Pendulum") + + axs[0,0].plot(self.theta) + axs[0,0].set_title('Angle [rad]') + axs[0,1].plot(self.v_ang) + axs[0,1].set_title('Angular velocity [rad/s]') + axs[1,0].plot(self.a_ang) + axs[1,0].set_title('Angular acceleration [rad/s^2]') + + fig, axs = plt.subplots(2, 2) + fig.suptitle("Cart") + + axs[0,0].plot(self.s_cart) + axs[0,0].set_title('Position [m]') + axs[0,1].plot(self.v_cart) + axs[0,1].set_title('Speed [m/s]') + axs[1,0].plot(self.a_cart) + axs[1,0].set_title('Acceleration [m/s^2]') + + plt.show() + + def pidControl(self, dt): + error = self.theta[self.index] + # (dt/1000) + # result = (self.kp * error) + (self.ki * sum(self.theta)) + (self.kd * self.v_ang[self.index]) #PID + result = (self.kp * error) + (((error - self.theta[self.index - 1]) / (dt + 1 / 1000)) * self.kd) #PD + self.a_cart[self.index] = result * 10 + \end{lstlisting} + \subsection{uiHelpers.py} + \begin{lstlisting}[language=Python] + import pygame + + # Constants + C_GRID_L_VALUE = 200 + C_GRID_D_VALUE = 100 + C_MPLOT_START = 700 + C_BLINK_TIME = 500 + + gridLight = pygame.Color(C_GRID_L_VALUE, C_GRID_L_VALUE, C_GRID_L_VALUE) + gridDark = pygame.Color(C_GRID_D_VALUE, C_GRID_D_VALUE, C_GRID_D_VALUE) + font_h = pygame.font.SysFont(None, 28) + font_m = pygame.font.SysFont(None, 16) + + # UI Class + class SimUI: + def __init__(self, screen, pole): + self.screen = screen + self.pole = pole + + self.metaPlotY = C_MPLOT_START + + self.blink = False + self.blinkTimer = 0 + + self.controlled = True + self.paused = False + + def meta(self, val, desc): + # Capture values to display status + if desc == "Control": + self.controlled = val + if desc == "Paused": + self.paused = val + + # Print em + self.screen.blit( + font_m.render(f"{desc} = {val}", True, "black"), (10, self.metaPlotY) + ) + self.metaPlotY -= 15 + + def grid(self, dist, Xoff=0, Yoff=0): + # Clear the screen + self.screen.fill("white") + + # Drawing offsets so the grid aligns with the pole + cXoff = self.pole.x % dist + cYoff = self.pole.y % dist + + # Draw the grid + for i in range(0, 1280, dist): + pygame.draw.line( + self.screen, + gridLight, + (i + Xoff + cXoff, 0), + (i + Xoff + cXoff, 720), + 1, + ) + pygame.draw.line( + self.screen, + gridLight, + (0, i + Yoff + cYoff), + (1280, i + Yoff + cYoff), + 1, + ) + + # Draw the center lines darker + pygame.draw.line( + self.screen, gridDark, (self.pole.x + Xoff, 0), (self.pole.x + Xoff, 720), 1 + ) + pygame.draw.line( + self.screen, + gridDark, + (0, self.pole.y + Yoff), + (1280, self.pole.y + Yoff), + 1, + ) + + def centeredText(self, font, text="", colour="black", y=0): + textObj = font.render(text, True, colour) + text_rect = textObj.get_rect(center=(1280 / 2, 720 / 2 - y)) + self.screen.blit(textObj, text_rect) + + def gameover(self, time): + font_g = pygame.font.Font("res\\pricedown.otf", 128) + self.centeredText(font_g, "wasted", "red", 150) + self.centeredText(font_m, f"You controlled the pendulum for {time / 1000} seconds", "black", 80) + + self.centeredText(font_m, "Press space to restart", "black", 60) + self.centeredText(font_m, "Press G to view nerd graphs", "black", 45) + + def update(self, dt): + # Credits + self.screen.blit( + font_h.render("Pendulum simulator 4000", True, "black"), (10, 10) + ) + self.screen.blit( + font_m.render("Arne van Iterson, 2023", True, "black"), (1150, 700) + ) + + # Reset meta printing + self.metaPlotY = C_MPLOT_START + + # Draw current status blinking + if self.blink: + if self.paused: + text = "Paused" + elif self.controlled: + text = "Auto mode" + else: + text = "Manual mode" + + textObj = font_h.render(text, True, "black") + text_rect = textObj.get_rect(right=1270, top=10) + self.screen.blit(textObj, text_rect) + + if self.blinkTimer > C_BLINK_TIME: + self.blink = not self.blink + self.blinkTimer = 0 + self.blinkTimer += dt + + \end{lstlisting} + \newpage + \section{Control source code} + \begin{lstlisting}[language=Python] + #!/usr/bin/env python3 + + # https://github.com/ev3dev/ev3dev-lang-python + # https://github.com/ev3dev/ev3dev-lang-python-demo#balanc3r + + from ev3dev2.motor import LargeMotor, OUTPUT_B, OUTPUT_C + from ev3dev2.sensor import INPUT_2 + from ev3dev2.sensor.lego import GyroSensor + from ev3dev2.sound import Sound + + import time + import logging + import os + + from datetime import datetime + + now = datetime.now() # current date and time + LOGNAME = os.path.join(os.getcwd(), "logs/", (now.strftime("%m-%d_%H:%M:%S") + ".csv")) + # print(LOGNAME) + logging.basicConfig(filename=LOGNAME, + filemode='a', + format='%(message)s', + level=logging.DEBUG) + + # Motors + l_motor = LargeMotor(OUTPUT_B) + r_motor = LargeMotor(OUTPUT_C) + + # Sensors + gyro = GyroSensor(INPUT_2) + + # Sound + sound = Sound() + sound.speak('3, 2, 1') + + def reset_gyro(): + gyro.calibrate() + gyro.mode = GyroSensor.MODE_GYRO_G_A + # if gyro.mode == gyro.MODE_GYRO_CAL: + # gyro.mode = gyro.MODE_GYRO_G_A + # gyro.mode = gyro.MODE_GYRO_CAL + # else: + # mode = gyro.mode + # gyro.mode = gyro.MODE_GYRO_CAL + # gyro.mode = mode + + # tires are 56mm diameter + def balance(target_angle=0, k_p=150, k_i=0, k_d=0.1): + integrated_error = 0 + t = time.time() + + # Calibrate gyro in current position + reset_gyro() + angle = gyro.angle + + # Stop if the robot has fallen over + while abs(angle) < 40: + # Keep time + dt = time.time() - t + t = time.time() + + angle, rate = gyro.angle_and_rate + error = angle - target_angle + integrated_error = integrated_error + error * dt + + pid = k_p * error + k_i * integrated_error + k_d * rate + + # print(t, angle, integrated_error, rate, pid) + logging.debug("%d;%d;%d", angle, rate, pid) + + if abs(pid) > 1050: + print("cap") + + speed = min(max(pid, -1050), 1050) + l_motor.run_forever(speed_sp=speed) + r_motor.run_forever(speed_sp=speed) + + l_motor.stop(stop_action="hold") + r_motor.stop(stop_action="hold") + + balance() + \end{lstlisting} + \end{document} % This is the end of the document diff --git a/doc/plot.ipynb b/doc/plot.ipynb new file mode 100644 index 0000000..30cff79 --- /dev/null +++ b/doc/plot.ipynb @@ -0,0 +1,68 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAHNCAYAAAAda3+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABj9UlEQVR4nO3de1xUdf4/8NcwMDN44aLcEUToguWFFhXxkpVstLqW39wiNS9kuhb0M2lLyQvZRaw1lzYpvprmtqthuupW+sU1ksxELdTySikopg6CyIAgt5nP7w+aoyMXucwNzuv5eMyj5jOfc877oPPxzed8LgohhAARERGRTDnYOgAiIiIiW2IyRERERLLGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLLGZIiILGLdunVQKBTSy9HREf7+/pg+fTouXLhgs7jeeustPProo/D29oZCocBrr71ms1iIyD442joAIurcXn/9dfTp0wdVVVXYv38/1q1bh7179+LYsWPQaDRWj2fhwoXw8fHBfffdh507d1r9+kRkf5gMEZFF/eEPf8CgQYMAAM8++yw8PDzw9ttv4/PPP8eTTz5p9Xjy8/MRFBSE4uJieHp6Wv36RGR/+JiMiKxq5MiRAIAzZ85IZQ888AAeeOCBBnWnT5+OoKAg6f3Zs2ehUCiwfPlyrFq1CiEhIVCr1Rg8eDC+//77Fl3/5vMREQHsGSIiKzt79iwAwN3dvc3n2LBhA8rLy/HnP/8ZCoUC77zzDh5//HHk5eXBycnJTJESkVwwGSIii9LpdCguLkZVVRUOHDiAJUuWQK1W449//GObz1lQUIBffvlFSqjuvvtuPPbYY9i5c2e7zktE8sRkiIgsKioqyuR9UFAQ/vWvf6FXr15tPmdMTIxJz5Lx0VteXl6bz0lE8sVkiIgsKjU1FXfddRd0Oh3Wrl2LPXv2QK1Wt+ucgYGBJu+NidHVq1fbdV4ikicmQ0RkUUOGDJFmk40fPx4jRozApEmTkJubi27dugEAFAoFhBANjtXr9Y2eU6lUNlre2DmIiG6Hs8mIyGqUSiWSk5Nx8eJFrFy5Uip3d3dHaWlpg/rnzp2zYnREJFdMhojIqh544AEMGTIEKSkpqKqqAgCEhITg1KlTKCoqkur9+OOP+O6772wVJhHJCB+TEZHVvfzyy3jiiSewbt06zJ49G8888wxWrFiB6OhozJgxA5cvX0ZaWhruvfdelJWVmfXa//znP3Hu3DlUVlYCAPbs2YM333wTADBlyhT07t3brNcjIvvHZIiIrO7xxx9HSEgIli9fjpkzZ6Jv37745JNPsHjxYiQkJOCee+7BP//5T2zYsAFZWVlmvfaaNWvwzTffSO93796N3bt3AwBGjBjBZIhIhhSCIw6JiIhIxjhmiIiIiGSNyRARERHJGpMhIiIikjUmQ0RERCRrTIaIiIhI1pgMERERkawxGSIiIiJZYzJEREREssZkiIiIiGSNyRARERHJGpMhIiIikjUmQ2QXFAoFXnvtNbOc6+zZs1AoFNJr8+bNZjnv7bz22mtQKBQmZW5ublIc8fHxVomDiIhah8kQmfjggw+gUCgQERFh61DabdasWfjnP/+JIUOG2CyGVatW4Z///KfNrk9ERLfHZIhMrF+/HkFBQTh48CBOnz5t63DaJTIyEk8//TQCAwNtFsOTTz6Jp59+2mbXJyKi22MyRJL8/Hzs27cPK1asgKenJ9avX2/rkKyioqLC1iEQEZENMRkiyfr16+Hu7o6xY8fiT3/6U6PJkHE8zvLly7Fq1SqEhIRArVZj8ODB+P777xvU37RpE+655x5oNBr069cPW7duxfTp0xEUFHTbeC5cuIBnnnkG3t7eUKvVuPfee7F27dp23aNxXM+JEycwadIkuLu7Y8SIEQCAn376CdOnT0dwcDA0Gg18fHzwzDPP4MqVKw3Os3fvXgwePBgajQYhISH43//933bFRUREtuNo6wDIfqxfvx6PP/44VCoVJk6ciA8//BDff/89Bg8e3KDuhg0bUF5ejj//+c9QKBR455138PjjjyMvLw9OTk4AgO3btyMmJgb9+/dHcnIyrl69ihkzZsDf3/+2sRQWFmLo0KHSwGNPT0/83//9H2bMmIGysjK8+OKL7brXJ554AnfeeSeWLl0KIQQAYNeuXcjLy0NsbCx8fHxw/PhxrFq1CsePH8f+/fulwdFHjx7Fww8/DE9PT7z22muoq6tDUlISvL292xUTERHZiCASQvzwww8CgNi1a5cQQgiDwSB69eol5syZY1IvPz9fABA9e/YUJSUlUvl//vMfAUB88cUXUln//v1Fr169RHl5uVSWlZUlAIjevXubnBeASEpKkt7PmDFD+Pr6iuLiYpN6Tz31lHB1dRWVlZVN3osxxo8//rjBZ0lJSQKAmDhxYoPPGjvnp59+KgCIPXv2SGXjx48XGo1GnDt3Tio7ceKEUCqVoqmvFAARFxfXZMxERGQ7fExGAOp7hby9vfHggw8CqJ/qHhMTg/T0dOj1+gb1Y2Ji4O7uLr0fOXIkACAvLw8AcPHiRRw9ehRTp05Ft27dpHqjRo1C//79m41FCIF///vfGDduHIQQKC4ull7R0dHQ6XQ4dOhQu+539uzZDcqcnZ2l/6+qqkJxcTGGDh0KANL19Ho9du7cifHjx5sMzO7bty+io6PbFRMREdkGkyGCXq9Heno6HnzwQeTn5+P06dM4ffo0IiIiUFhYiMzMzAbH3DpDy5gYXb16FQBw7tw5AMAdd9zR4NjGym5WVFSE0tJSrFq1Cp6eniav2NhYAMDly5dbf6M36dOnT4OykpISzJkzB97e3nB2doanp6dUT6fTSbFdv34dd955Z4Pj77777nbFREREtsExQ4Svv/4aly5dQnp6OtLT0xt8vn79ejz88MMmZUqlstFzid/G37SHwWAAADz99NOYNm1ao3UGDBjQrmvc3Atk9OSTT2Lfvn14+eWXERYWhm7dusFgMOCRRx6RYiIios6HyRBh/fr18PLyQmpqaoPPtmzZgq1btyItLa3RBKIpvXv3BoBG1yq63fpFnp6e6N69O/R6PaKiolp8zfa4evUqMjMzsWTJEixevFgq/+WXXxrE5uzs3KAcAHJzcy0eJxERmR8fk8nc9evXsWXLFvzxj3/En/70pwav+Ph4lJeX4/PPP2/Vef38/NCvXz988sknuHbtmlT+zTff4OjRo80eq1QqMWHCBPz73//GsWPHGnxeVFTUqlhawtjTdWvPVkpKSoN60dHR2LZtGwoKCqTykydPYufOnWaPi4iILI89QzL3+eefo7y8HI8++mijnw8dOlRagDEmJqZV5166dCkee+wxDB8+HLGxsbh69SpWrlyJfv36mSRIjVm2bBl2796NiIgIzJw5E/fccw9KSkpw6NAhfPXVVygpKWlVLLfj4uKC+++/H++88w5qa2vh7++P//73v8jPz29Qd8mSJcjIyMDIkSPx/PPPo66uDu+//z7uvfde/PTTT2aNi4iILI89QzK3fv16aDQa/P73v2/0cwcHB4wdOxYZGRmNLj7YnHHjxuHTTz9FTU0N5s+fjy1btmDdunW4++67odFomj3W29sbBw8eRGxsLLZs2YL4+Hi89957KCkpwdtvv92qOFpqw4YNiI6ORmpqKhITE+Hk5IT/+7//a1BvwIAB2LlzJzw9PbF48WKsXbsWS5Yswf/8z/9YJC4iIrIshTDHiFeiVggLC4Onpyd27dplkfOfPXsWffr0wfvvv4+nnnoKLi4uUKlUFrnW7ZSUlMBgMMDT0xNxcXFYuXKlTeIgIqKmsWeILKa2thZ1dXUmZVlZWfjxxx/xwAMPWPz6L7zwAjw9PVs93smcgoOD4enpabPrExHR7bFniCzm7NmziIqKwtNPPw0/Pz+cOnUKaWlpcHV1xbFjx9CzZ0+LXLeqqgp79+6V3g8YMABeXl4WudbtfPPNN6itrQUABAQEcC0iIiI7xGSILEan02HWrFn47rvvUFRUhK5du2L06NFYtmwZQkJCbB0edWB79uzBX//6V+Tk5ODSpUvYunUrxo8f3+wxWVlZSEhIwPHjxxEQEICFCxdi+vTpVomXiOwbZ5ORxbi6umLjxo22DoM6oYqKCgwcOBDPPPMMHn/88dvWz8/Px9ixYzF79mysX78emZmZePbZZ+Hr68ttVIiIPUNE1LEpFIrb9gzNmzcP27dvN1m36qmnnkJpaSkyMjKsECUR2TMOoCaiTi87O7vBaubR0dHIzs62UUREZE9k+ZjMYDDg4sWL6N69OxQKha3DIZIlIQTKy8vh5+cHBwfL/l6m1Wrh7e1tUubt7Y2ysjJcv369ya1mqqurUV1dLb03GAwoKSlBz5492XYQ2YCl2g1ZJkMXL15EQECArcMgIgDnz59Hr169bB1Go5KTk7FkyRJbh0FEtzB3uyHLZKh79+4A6n+YLi4uNo6GSJ7KysoQEBAgfR8tycfHB4WFhSZlhYWFcHFxaXYD4sTERCQkJEjvdTodAgMD2XYQ2Yil2g1ZJkPG7m0XFxc2aEQ2Zo3HTZGRkdixY4dJ2a5duxAZGdnscWq1Gmq1ukE52w4i2zJ3u8EB1ETU4Vy7dg1HjhzBkSNHANRPnT9y5AgKCgoA1PfoTJ06Vao/e/Zs5OXl4ZVXXsGpU6fwwQcf4LPPPsPcuXNtET4R2RkmQ0TU4fzwww+47777cN999wEAEhIScN9992Hx4sUAgEuXLkmJEQD06dMH27dvx65duzBw4EC8++67+Oijj7jGEBEBkOk6Q2VlZXB1dYVOp2NXN5GNdMTvYUeMmagzsdR3kD1DRHauqlaP+f/+CTuPa20dChFRp8RkiMjOZeUWIf3783g745StQyEi6pSskgylpqYiKCgIGo0GEREROHjwYLP1N23ahNDQUGg0GvTv37/BLJCbzZ49GwqFAikpKWaOmsg+nCm6BgAouFKJWr3BxtEQEXU+Fk+GNm7ciISEBCQlJeHQoUMYOHAgoqOjcfny5Ubr79u3DxMnTsSMGTNw+PBhjB8/HuPHjzfZU8ho69at2L9/P/z8/Cx9G0Q2k1dUAQCoMwgUlFTaOBoios7H4snQihUrMHPmTMTGxuKee+5BWloaunTpgrVr1zZa/7333sMjjzyCl19+GX379sUbb7yB3/3ud1i5cqVJvQsXLuCFF17A+vXr4eTkZOnbILKZvOJrN/7/t8SIiIjMx6LJUE1NDXJyckw2SHRwcEBUVFSTGyS2ZENFg8GAKVOm4OWXX8a9995rmeCJ7IAQwiQByiu61kxtIiJqC4uuQF1cXAy9Xt/oBomnTjU+GLSpDRW12hszad5++204Ojri//2//9eiOG7dbLGsrKylt0BkU1cra6G7Xiu9Z88QEZH5dbjZZDk5OXjvvfewbt26Fi/HnZycDFdXV+nFTVqpo7i1J+jmR2ZERGQeFk2GPDw8oFQqG90g0cfHp9FjmtpQ0Vj/22+/xeXLlxEYGAhHR0c4Ojri3LlzeOmllxAUFNToORMTE6HT6aTX+fPn239zRFZg7Anq2VUFAMgvZs8QEZG5WTQZUqlUCA8PR2ZmplRmMBiQmZnZ5AaJkZGRJvUB0w0Vp0yZgp9++knal+jIkSPw8/PDyy+/jJ07dzZ6TrVaLW2syA0WqSM581tP0AN3ewEAiq/VmDw2IyKi9rP4rvUJCQmYNm0aBg0ahCFDhiAlJQUVFRWIjY0FAEydOhX+/v5ITk4GAMyZMwejRo3Cu+++i7FjxyI9PR0//PADVq1aBQDo2bMnevbsaXINJycn+Pj44O6777b07RBZlbFnqL+/C779RY3L5dXIK7qG+wLdbRwZEVHnYfFkKCYmBkVFRVi8eDG0Wi3CwsKQkZEhDZIuKCiAg8ONDqphw4Zhw4YNWLhwIV599VXceeed2LZtG/r162fpUInsjnHMULBnNwR7dv0tGapgMkREZEYWT4YAID4+HvHx8Y1+lpWV1aDsiSeewBNPPNHi8589e7aNkRHZrzq9QVpkMdizK4I9u2F/XgkHURMRmZlVkiEiMrVg61GsP1DQoroaJwf4uToj2KMrAE6vJyIytw43tZ6oozMYBP596NcW1x/d1xsODgqEeHYDwGSIiMjc2DNEZGWXyqpQVWuAo4MC+xIfgrKZ9bIUCgXcu9RvNxPsWd8zlH+lAnqDgNKhZetsERFR85gMEVlZ/m89O4E9u8Cru6bFx/Vy7wInpQI1dQZcLL2OgB5dLBUiEZGs8DEZkZUZB0AHe3Rr1XFKBwV69/xt3BAXXyQiMhsmQ0RWZhzzE/LbY6/WuDGImjPKiIjMhckQkZWdkdYOakMyxEHURERmx2SIyMqMiYwxsWkNYwLFtYaIiMyHyRCRFVXV6nFRdx0A0Mej9T1Dxkdr7BkiIjIfJkNEVpRfXAEhABeNo7QTfWsYB11f0lWhsqbO3OEREckSkyEiK8ovvvGITNHM+kJNce+qktYdYu8QkJqaiqCgIGg0GkRERODgwYPN1k9JScHdd98NZ2dnBAQEYO7cuaiqqrJStERkr5gMEVlRXjsGTxtJg6hlPr1+48aNSEhIQFJSEg4dOoSBAwciOjoaly9fbrT+hg0bMH/+fCQlJeHkyZNYs2YNNm7ciFdffdXKkRORvWEyRGQF12v02He6GAfPXgUAaWuNtuD0+norVqzAzJkzERsbi3vuuQdpaWno0qUL1q5d22j9ffv2Yfjw4Zg0aRKCgoLw8MMPY+LEibftTSKizo/JEJEVzN/yEyZ9dAB7fi4CcCOhaYs+HESNmpoa5OTkICoqSipzcHBAVFQUsrOzGz1m2LBhyMnJkZKfvLw87NixA2PGjGnyOtXV1SgrKzN5EVHnw+04iKwg51x9j1CwZ1cEe3TF/Xd5tvlcxkHU+TJ+TFZcXAy9Xg9vb2+Tcm9vb5w6darRYyZNmoTi4mKMGDECQgjU1dVh9uzZzT4mS05OxpIlS8waOxHZH/YMEVlYVa0eF0rrp9N/9udIfDRtMLqq2/57yI3p9dcghDBLjHKQlZWFpUuX4oMPPsChQ4ewZcsWbN++HW+88UaTxyQmJkKn00mv8+fPWzFiIrIW9gwRWdjZK+2bTn+rwJ5d4KAAKmr0uFxeDW+Xlm/22ll4eHhAqVSisLDQpLywsBA+Pj6NHrNo0SJMmTIFzz77LACgf//+qKiowKxZs7BgwQI4ODT83VCtVkOtVpv/BojIrrBniMjCbl5xui3T6W+ldlRKO9afkekgapVKhfDwcGRmZkplBoMBmZmZiIyMbPSYysrKBgmPUqkEAPawEckckyEiC5Om07dj0PStbswok++4oYSEBKxevRr/+Mc/cPLkSTz33HOoqKhAbGwsAGDq1KlITEyU6o8bNw4ffvgh0tPTkZ+fj127dmHRokUYN26clBQRkTxZJRlq7cJomzZtQmhoKDQaDfr3748dO3ZIn9XW1mLevHno378/unbtCj8/P0ydOhUXL1609G0QtcmNniEzJkPcsBUxMTFYvnw5Fi9ejLCwMBw5cgQZGRnSoOqCggJcunRJqr9w4UK89NJLWLhwIe655x7MmDED0dHR+N///V9b3QIR2QmLJ0OtXRht3759mDhxImbMmIHDhw9j/PjxGD9+PI4dOwagvqv70KFDWLRokTQIMjc3F48++qilb4WoTc4Ut31j1qZww9Z68fHxOHfuHKqrq3HgwAFERERIn2VlZWHdunXSe0dHRyQlJeH06dO4fv06CgoKkJqaCjc3N+sHTkR2RSEs/LA8IiICgwcPxsqVKwHUP9cPCAjACy+8gPnz5zeoHxMTg4qKCnz55ZdS2dChQxEWFoa0tLRGr/H9999jyJAhOHfuHAIDA28bU1lZGVxdXaHT6eDi4tLGOyO6PSEEBiz5L8qr6pDx4kiE+pjn71v2mSuYuHo/Ant0wZ5XHjTLOa2tI34PO2LMRJ2Jpb6DFu0ZasvCaNnZ2Sb1ASA6OrrJ+gCg0+mgUCia/A2PC6eRrRRfq0F5VR0UCiCop/kekxmn1/96tRLVdXqznZeISI4smgw1tzCaVqtt9BitVtuq+lVVVZg3bx4mTpzYZJaYnJwMV1dX6RUQENCGuyFqPePCiP5uztA4mW+Qrmd3NbqqlDAIoOBKpdnOS0QkRx16NlltbS2efPJJCCHw4YcfNlmPC6eRrRhnkvUx40wyAFAoFNIYpDMyHkRNRGQOFl10sS0Lo/n4+LSovjEROnfuHL7++utmnx1y4bTOY3POr8jVtv4xZ3eNE6YPD4KLxskCUTXNuLN8ezZmbUqwZ1ccvaCT/SBqIqL2smgydPPCaOPHjwdwY2G0+Pj4Ro+JjIxEZmYmXnzxRals165dJgupGROhX375Bbt370bPnj0teRtkJ/KLK/CXTT+2+Xi1owP+PCrEjBHdnrTGkBmn1RsZ9yiT8/R6IiJzsPh2HAkJCZg2bRoGDRqEIUOGICUlpcHCaP7+/khOTgYAzJkzB6NGjcK7776LsWPHIj09HT/88ANWrVoFoD4R+tOf/oRDhw7hyy+/hF6vl8YT9ejRAypV+7c7IPtk7BHyddXg0YF+LT7up191yM67glxtuaVCa5K0xpCHZXqG6q/BniEiovaweDIUExODoqIiLF68GFqtFmFhYQ0WRrt5ifxhw4Zhw4YNWLhwIV599VXceeed2LZtG/r16wcAuHDhAj7//HMAQFhYmMm1du/ejQceeMDSt0Q2YhwbE9GnBxLH9G3xcTuOXkJ23hVpvR9rqdUbUFBSP7jZIj1D0lpD7BkiImoPq2zUGh8f3+RjsaysrAZlTzzxBJ544olG6wcFBXEfIZm6eY+v1gi+ZZd3c+wP1hIFJZWoMwg4OynhY4HNVI2Dsksra1FSUYMeZtgElohIjjr0bDKSF+NA4db2sgT17AqFAiivqkPxtRpLhNYoY/LWx6MrHBzMn4B1UTnCz1Xz27X4qIyIqK2YDFGHYVyzp7XjbzROSvi7OQOwbtJgycHTRtIeZXxURkTUZkyGqEMoqahBaWUtgLat2WOLpOHG4GnLJUN9uHs9EVG7MRmiDsHYy+Lv5gxnVetXcjYmJPlWTIbyLbBB6604o4yIqP2YDFGHcGPwdNt6WUJskDS0dYxTa/AxGRFR+zEZog7hTHH7trXoY+UFCnXXa6XB2ubeiuNmxh6vc1cqUKc3WOw6RESdGZMh6hDaO/7G2DtTUFKJWiskDcYeKK/uanS34BYg/m7OUDs6oFYv8OvV6xa7DhFRZ8ZkiDqEGzOz2jb+xsdFA2cnJeoMQloI0ZLa+1ivpRwcFDcGUXOPMiKiNrHKootELfXRt3nYnPNrg/Ibg5HbllwYk4YTl8qQV1RhkY1Tb3ZjvJBlr1N/ja44pS3Hq1uOIaDHGSQ/PgB3eFn+ukREnQV7hshuCCHw3le/4JS2vMHLIOr3JPNzdW7z+a0588oa0+qNwnv3AABoy6rw/dmrjSaTRETUNPYMkd0oKq9GeXUdHBTAx7FDoLxl24xQ3+7tWslZmnllhUHUxmtYugcKAKYPC8J9gW74/MhFrNt3ltPsiYhaickQ2Q3j9PBe7l0w6i5Ps5/fOL3e0msNGQwC+VesM2YIAJQOCvwu0B3lVXX1yRCn2RMRtQofk5HdsPSgY2sNNL5Qeh01dQY4KRXo5d7Fote62c3T7PUGbmZMRNRSTIbIbkgzxlq591hLGZOh4ms10F2vtcg1gBs9XL17doXSAhu0NsV0mr3lZ8wREXUWTIbIbuS1c8bY7XTXOMGru7r+WhYcV3MjqbP8I7KbmUyzl8leZampqQgKCoJGo0FERAQOHjzYbP3S0lLExcXB19cXarUad911F3bs2GGlaInIXjEZIrthjSTixowyyyULNx73WX96uzEZOiODQdQbN25EQkICkpKScOjQIQwcOBDR0dG4fPlyo/Vramrw+9//HmfPnsXmzZuRm5uL1atXw9/f38qRE5G94QBqsgs1dQac/20FZctubNoN+/NKLDpuyBp7kjVFSvZkMIh6xYoVmDlzJmJjYwEAaWlp2L59O9auXYv58+c3qL927VqUlJRg3759cHKqXxU8KCjImiETkZ1izxDZhYKSSugNAl1VSni7qC12HWvsXn9jWr0NkiFpD7bO3TNUU1ODnJwcREVFSWUODg6IiopCdnZ2o8d8/vnniIyMRFxcHLy9vdGvXz8sXboUer2+yetUV1ejrKzM5EVEnY9VkqHWPtfftGkTQkNDodFo0L9//wbP9IUQWLx4MXx9feHs7IyoqCj88ssvlrwFsjDjP959PLtCobDcoOMQC681VFlTh0u6KgCWGwjeHGs8BrQHxcXF0Ov18Pb2Nin39vaGVqtt9Ji8vDxs3rwZer0eO3bswKJFi/Duu+/izTffbPI6ycnJcHV1lV4BAQFmvQ8isg8WT4Za+1x/3759mDhxImbMmIHDhw9j/PjxGD9+PI4dOybVeeedd/D3v/8daWlpOHDgALp27Yro6GhUVVVZ+nbIQqTB0xZOIIJvWmvIYIHp58YeJ/cuTnDvqjL7+W/H+Ijxcnk1yqssN2OuIzIYDPDy8sKqVasQHh6OmJgYLFiwAGlpaU0ek5iYCJ1OJ73Onz9vxYiJyFosPmaotc/133vvPTzyyCN4+eWXAQBvvPEGdu3ahZUrVyItLQ1CCKSkpGDhwoV47LHHAACffPIJvL29sW3bNjz11FNmi113vRY1dZbf4byzUzk6wNW56Z3bSypqcOpS/eMHS4+z6eXeBSqlA6rrDDh+sQx3eneDxklptvPbcvA0ALg6O8GjmwrF12qQX1yBAb3cbBKHpXl4eECpVKKwsNCkvLCwED4+Po0e4+vrCycnJyiVN/68+/btC61Wi5qaGqhUDZNXtVoNtdpyj22JyD5YNBkyPtdPTEyUym73XD87OxsJCQkmZdHR0di2bRsAID8/H1qt1mSsgKurKyIiIpCdnd1oMlRdXY3q6mrpfUuf+7/02RF8dbLxHixqneTH+2PikMAG5W9tP4HV3+ZL7y2dRCgdFOjdswt+uXwN41buRTe1I3bOvR/+bm3f8+xm1tyTrCnBHt1QfK0EeUWdNxlSqVQIDw9HZmYmxo8fD6C+5yczMxPx8fGNHjN8+HBs2LABBoMBDg71neI///wzfH19G02EiEg+LPqYrC3P9bVabbP1jf9tzTnb/txfAYUCfLXjZZTZRFJ5c7Lp7+aMYSE9W/hn03aP/64XnJT1wV2rrsOBvCtmO3dhef2jWj8zJVdtcWOl7c49bighIQGrV6/GP/7xD5w8eRLPPfccKioqpF7oqVOnmvwi9txzz6GkpARz5szBzz//jO3bt2Pp0qWIi4uz1S0QkZ2QxdT6xMREk96msrKyFiVEH00bZMmwZGHvL8V4es2BRqey19QZUFBSv1Ly/sTR8HHVWCWm5x4IwXMPhODVrUex4UCBWQcbl1bWAAB62GC8kNGNQdSde0ZZTEwMioqKsHjxYmi1WoSFhSEjI0P6RamgoEDqAQKAgIAA7Ny5E3PnzsWAAQPg7++POXPmYN68eba6BSKyExZNhtryXN/Hx6fZ+sb/FhYWwtfX16ROWFhYo+fkc3/bMf7DXHClErV6A5yUN/5xMk6n72Lh6fRNxmaBvcpKKuqTIVsMnjYKtvCMOXsSHx/f5GOxrKysBmWRkZHYv3+/haMioo7Goo/Jbn6ub2R8rh8ZGdnoMZGRkSb1AWDXrl1S/T59+sDHx8ekTllZGQ4cONDkOcl2fFw00Dg5oM4gcL7EdL8saTq9h2Wn0zfFEtPQSyvrZ3C5d2l6wLilWXrGHBFRZ2PxqfWtfa4/Z84cZGRk4N1338WpU6fw2muv4YcffpB++1MoFHjxxRfx5ptv4vPPP8fRo0cxdepU+Pn5SQMpyX7U75fVeE/Fjb3IbDPzyjiN35xJg9Qz1MV2PUOBPbrA0UGB67V6aMu43AQR0e1YfMxQa5/rDxs2DBs2bMDChQvx6quv4s4778S2bdvQr18/qc4rr7yCiooKzJo1C6WlpRgxYgQyMjKg0VhnzAm1TrBnV5y8VPbb46gbA99ttaGpUS93ZzgpFaiuM+BC6XUE9OjSrvMJIW70DNnwMZmT0gGBPbogr7gCeUUVNh3MTUTUEVhlAHVrn+s/8cQTeOKJJ5o8n0KhwOuvv47XX3/dXCGSBYU0sQVGvoV3qb8dR6UDevfsitOXryGvuKLdyVBljR41+vp1qWz5mAyo/5nmFVcgr/gaRtzpYdNYiIjsHfcmI4szPgY7c+tjMmkPL9s8JgNuGkRthplXxkdkakcHOJtxIce2kNMgaiKi9mIyRBbX2EBlXWUtrvyWPPSx5QKFZkwabgyeVtlkQPjNjD/TM518ej0RkTkwGSKLM/7DXHytGmW/7Zd15rfp7D4uGnRV2265KylRM8P0+pJK20+rN7rR48WeISKi22EyRBbXXeMEz+716wgZ/3E2/teWvULAjaQh3yw9Q8aZZLYdLwTc6PG6qLuOqlq9jaMhIrJvTIbIKoxJx5ZDv2LLoV+x60T91im2GjxtdCNpqEJlTV27zmUPCy4aeXRTobvGEUIAZ6+wd4iIqDlMhsgq7vCqTzo+yT6HhM9+xM7j9auM23LwNFC/bYbbbz05t852a62rdrDgopFCoeAgaiKiFpLF3mRke7HDg3DlWg0qb3pk497FCePv87dhVPWCPbriUEEp8ooqcK+fa5vPc9UOFly8WYhHV/x4vrTT71FGRNReTIbIKu7w6o60KeG2DqNRwZ7dpGSoPa5W2lcyZIntRoiIOiM+JiPZM9eMshurT9v+MRlw0/pO7Xz8R0TU2TEZItkLbmLvtNayh33JbiZt2Fp0DUJww1YioqYwGSLZu3mX9/YkDaV29pgsqGdXKBRAWVWdtMAlERE1xGSIZK93zy5wUADXqutQVF7d5vMYF13sYQdT6wFA46SEn2v9Jq0cN0RE1DQmQyR7akclernXb9J66/5pLXW9Ro+q2vpNWt3sYGq90Y1B1JxRRkTUFCZDRGj/IGrjTDInpQLdbLi9yK2M6zjlcRA1EVGTmAwRof2DqI3JkJsdbNJ6M06vJyK6PSZDRGj/46SrFfaz+vTNpCTPDBvREhF1VkyGiGA6o6wt7G3BRSPjfRVcqUSt3mDjaIiI7JP9DG4gsiHj2JqCkkrM3XhEKh9xhwcmhPe67fH2Nq3eyMdFA2cnJa7X6nG+pFJaiNFSPsg6jbyiCjw1OACDgnpY9FpEROZi0Z6hkpISTJ48GS4uLnBzc8OMGTNw7Vrz3fVVVVWIi4tDz5490a1bN0yYMAGFhYXS5z/++CMmTpyIgIAAODs7o2/fvnjvvfcseRskA17d1fDopoJBAFsPX5Be8/79E67X6G97/OXfpuT36GZfyZCDgwJ9PKw3bijz5GVszvkVl3RVFr8WAKSmpiIoKAgajQYRERE4ePBgi45LT0+HQqHA+PHjLRsgEXUIFk2GJk+ejOPHj2PXrl348ssvsWfPHsyaNavZY+bOnYsvvvgCmzZtwjfffIOLFy/i8ccflz7PycmBl5cX/vWvf+H48eNYsGABEhMTsXLlSkveCnVyCoUC62KHYOHYvtKru9oRdQaBs1dun0QYZ2v16dnV0qG2mrm2G2kJ45gr4zUtaePGjUhISEBSUhIOHTqEgQMHIjo6GpcvX272uLNnz+Ivf/kLRo4cafEYiahjsNhjspMnTyIjIwPff/89Bg0aBAB4//33MWbMGCxfvhx+fn4NjtHpdFizZg02bNiAhx56CADw8ccfo2/fvti/fz+GDh2KZ555xuSY4OBgZGdnY8uWLYiPj7fU7ZAM9PN3RT//G7vWbz96CYd/28C1r69Ls8cae12skQS0VrBH+8ZDtdTVihpc/W1/NmNvlCWtWLECM2fORGxsLAAgLS0N27dvx9q1azF//vxGj9Hr9Zg8eTKWLFmCb7/9FqWlpRaPk4jsn8V6hrKzs+Hm5iYlQgAQFRUFBwcHHDhwoNFjcnJyUFtbi6ioKKksNDQUgYGByM7ObvJaOp0OPXo0PT6huroaZWVlJi+i27kx3b75HhWDQSC/2NgjYtkxOW0hbdhq4cdkxt4xX1cNuqgsOxyxpqYGOTk5Jm2Fg4MDoqKimm0rXn/9dXh5eWHGjBktug7bDiJ5sFgypNVq4eXlZVLm6OiIHj16QKvVNnmMSqWCm5ubSbm3t3eTx+zbtw8bN25s9vFbcnIyXF1dpVdAQEDrboZk6cbjpeaTiEtlVaiqNcDRQYEAd2drhNYq1lpryJqPyIqLi6HX6+Ht7W1S3lxbsXfvXqxZswarV69u8XXYdhDJQ6uTofnz50OhUDT7OnXqlCVibeDYsWN47LHHkJSUhIcffrjJeomJidDpdNLr/PnzVomPOraQFq49ZPw8sGcXOCrtb7UK4yOr4mvVKKuqtdh1jEmjsUfNnpSXl2PKlClYvXo1PDw8Wnwc2w4ieWh1X/ZLL72E6dOnN1snODgYPj4+DQYy1tXVoaSkBD4+Po0e5+Pjg5qaGpSWlpr0DhUWFjY45sSJExg9ejRmzZqFhQsXNhuPWq2GWq1utg7RrfrctCq1EKLJlaWl8UJ2mAQAQHeNE7y6q3G5vBp5RRUIC3CzyHWs2TPk4eEBpVJpMtMUaLytAIAzZ87g7NmzGDdunFRmMNSvu+To6Ijc3FyEhIQ0OI5tB5E8tDoZ8vT0hKen523rRUZGorS0FDk5OQgPDwcAfP311zAYDIiIiGj0mPDwcDg5OSEzMxMTJkwAAOTm5qKgoACRkZFSvePHj+Ohhx7CtGnT8NZbb7X2FohapHfPLlAogPLqOhRdq4ZXd02j9YwDk0PscPC0UbBn19+SoWsWTIaMg8gtnxSqVCqEh4cjMzNTmh5vMBiQmZnZ6ESK0NBQHD161KRs4cKFKC8vx3vvvcfHX0QyZ7FRjn379sUjjzyCmTNnIi0tDbW1tYiPj8dTTz0lzSS7cOECRo8ejU8++QRDhgyBq6srZsyYgYSEBPTo0QMuLi544YUXEBkZiaFDhwKofzT20EMPITo6GgkJCdL4AKVS2aIkjailNE5K9HJ3xvmS68gvqmgyGTpjxR6Rtgr27Ib9eSUWGzekNwicu1JZfy0rzCQDgISEBEybNg2DBg3CkCFDkJKSgoqKCml22dSpU+Hv74/k5GRoNBr069fP5Hhj7/Ot5UQkPxad8rF+/XrEx8dj9OjRcHBwwIQJE/D3v/9d+ry2tha5ubmorKyUyv72t79JdaurqxEdHY0PPvhA+nzz5s0oKirCv/71L/zrX/+Synv37o2zZ89a8nZIhoI9uuF8yXXkFVcgIrhno3Ws2SPSVsYExVJrDf16tRI1egPUjg7wd7POIPKYmBgUFRVh8eLF0Gq1CAsLQ0ZGhjSouqCgAA4O9jeGi4jsj0IIIWwdhLWVlZXB1dUVOp0OLi7Nrx9D8rbki+P4+LuzmDmyDxaMvafB51W1evRdnAEhgJyFUejZzT7Hl+w+dRmx675HqE93ZLx4v12cvyN+DztizESdiaW+g/y1iagZxt6eph4v5RdXQAjA1dkJPbra11YcN+tz08KLBoP5f/8xPiq0xmKLRETmxo1aiZphfLz0zc9FGPTmrgaf19TVz0jq49G1ydlm9qCXuzOclApU1xlwofQ6Anp0Mev5jYPI7XncFBFRU5gMETWjn58rXDSOKKuqQ/G1mibrDb+j8fFE9sJR6YDePbvi9OVryC+uMHsyZO/LCxARNYfJEFEzXLs44dt5D+GS7nqTdZyUDlabQdUewR71yVBe0TXcf5d5Z17mFdv/jDoioqYwGSK6DVdnJ7g6O9k6jHarH/9UeNvtRVrrWnUdCsuqb7oGEVHHwgHURDJhqT3K8n87n0c3VadIGolIfpgMEclES/daay3pERnHCxFRB8VkiEgmjMnKRV0VKmvqzHbeM0WcSUZEHRuTISKZcO+qgnuX+sdY+WYcN8Rp9UTU0TEZIpKR2y0i2RbSbvV8TEZEHRSTISIZuXklanMQQkjn6sOeISLqoJgMEclIsJkHUWvLqlBZo4ejgwKBZl7IkYjIWpgMEcmI8VGWudYaMj5uC+zRBU5KNidE1DGx9SKSkZCb1hoSon0btp67UoEDeVcAcPA0EXVsXIGaSEYCe3aBg6J+1eii8mp4uWjadJ5PDxYgcctR6T1Xniaijow9Q0QyonZUSpu0nmnHjLJ9Z+p7hLqrHRHs2RWPDvQzS3xERLbAniEimQn26IpzVyqRV3wNkSE923QO4wDsd58ciIfv9TFneEREVseeISKZae9aQzdPp+fjMSLqDCyaDJWUlGDy5MlwcXGBm5sbZsyYgWvXmp/SW1VVhbi4OPTs2RPdunXDhAkTUFhY2GjdK1euoFevXlAoFCgtLbXAHRB1Pu2dXm+cTq/kdHoi6iQsmgxNnjwZx48fx65du/Dll19iz549mDVrVrPHzJ07F1988QU2bdqEb775BhcvXsTjjz/eaN0ZM2ZgwIABlgidqNMyLrzY1un1N0+nVzmyc5mIOj6LtWQnT55ERkYGPvroI0RERGDEiBF4//33kZ6ejosXLzZ6jE6nw5o1a7BixQo89NBDCA8Px8cff4x9+/Zh//79JnU//PBDlJaW4i9/+YulboGoUwr57dHW+ZJK1NQZWn28MYkK9uB0eiLqHCyWDGVnZ8PNzQ2DBg2SyqKiouDg4IADBw40ekxOTg5qa2sRFRUllYWGhiIwMBDZ2dlS2YkTJ/D666/jk08+gYPD7W+huroaZWVlJi8iufLqrkZXlRIGARSUtL53SNqLjGsLEVEnYbFkSKvVwsvLy6TM0dERPXr0gFarbfIYlUoFNzc3k3Jvb2/pmOrqakycOBF//etfERgY2KJYkpOT4erqKr0CAgJaf0NEnYRCoZAGPrdler3xMRkHTxNRZ9HqZGj+/PlQKBTNvk6dOmWJWAEAiYmJ6Nu3L55++ulWHaPT6aTX+fPnLRYfUUcQfNNK1K2VV2zcpd72PUOpqakICgqCRqNBREQEDh482GTd1atXY+TIkXB3d4e7uzuioqKarU9E8tHqdYZeeuklTJ8+vdk6wcHB8PHxweXLl03K6+rqUFJSAh+fxtcl8fHxQU1NDUpLS016hwoLC6Vjvv76axw9ehSbN28GAGlLAQ8PDyxYsABLlixpcF61Wg21Wt3SWyTq9KQ9ylo5o6yqVo9fr16vP4eNe4Y2btyIhIQEpKWlISIiAikpKYiOjkZubm6DXmkAyMrKwsSJEzFs2DBoNBq8/fbbePjhh3H8+HH4+/vb4A6IyF60Ohny9PSEp6fnbetFRkaitLQUOTk5CA8PB1CfyBgMBkRERDR6THh4OJycnJCZmYkJEyYAAHJzc1FQUIDIyEgAwL///W9cv35dOub777/HM888g2+//RYhISGtvR0iWZJ6hlo5o+zclUoIAXTXOMKjm8oSobXYihUrMHPmTMTGxgIA0tLSsH37dqxduxbz589vUH/9+vUm7z/66CP8+9//RmZmJqZOnWqVmInIPllsBeq+ffvikUcewcyZM5GWloba2lrEx8fjqaeegp9f/dL9Fy5cwOjRo/HJJ59gyJAhcHV1xYwZM5CQkIAePXrAxcUFL7zwAiIjIzF06FAAaJDwFBcXS9e7dawRETWurWsN3Rg83Q0KhcLscbVUTU0NcnJykJiYKJU5ODggKirKZLJFcyorK1FbW4sePXo0Wae6uhrV1dXSe06+IOqcLLpIyPr16xEaGorRo0djzJgxGDFiBFatWiV9Xltbi9zcXFRWVkplf/vb3/DHP/4REyZMwP333w8fHx9s2bLFkmESyY5xraGrlbW4WlHT4uOMPUkhNh4vVFxcDL1eD29vb5Pymydb3M68efPg5+dnMnv1Vpx8QSQPFt2brEePHtiwYUOTnwcFBUljfow0Gg1SU1ORmpraoms88MADDc5BRM3ronKEn6sGF3VVyCu+hvCuTfeO3OxMJ5lWv2zZMqSnpyMrKwsajabJeomJiUhISJDel5WVMSEi6oS4USuRTPXx7IqLuiqcKapAeO+WJUPGPcn6eNh28LSHhweUSmWDrXpunmzRlOXLl2PZsmX46quvbruCPSdfEMkD19InkinjjLL8Fg6iFkLctMaQbXuGVCoVwsPDkZmZKZUZDAZkZmZKky0a88477+CNN95ARkaGyYKwRCRv7BkikqnWDqIuqaiB7notFIobY45sKSEhAdOmTcOgQYMwZMgQpKSkoKKiQppdNnXqVPj7+yM5ORkA8Pbbb2Px4sXYsGEDgoKCpLFF3bp1Q7duXECSSM6YDBHJlHGdoJYuvGgcPO3n6gyNk9JicbVUTEwMioqKsHjxYmi1WoSFhSEjI0MaVF1QUGCyXc+HH36Impoa/OlPfzI5T1JSEl577TVrhk5EdobJEJFMGVeQPnelEnqDgNKh+any9rgnWXx8POLj4xv9LCsry+T92bNnLR8QEXVIHDNEJFP+bs5QOzqgRm/Ar1crb1vf2IMUwj3JiKiTYTJEJFMODgpp7E9LHpWdsZPB00RE5sZkiEjGjInNmRYMor6xQSt7hoioc2EyRCRj0oatt5leX6s3oOBK/aM09gwRUWfDZIhIxm48Jmu+Z+h8SSXqDAIaJwf4uDS9YjMRUUfEZIhIxoy9PD+e12HKmgNYf+Bco/VuXnna4TazzoiIOhomQ0Qydqd3d3RRKXG9Vo9vfynGki9OQG9ouNefvaw8TURkCUyGiGSsm9oRW58fjpSYsPpp9nUGXLh6vUE94+BpW+9WT0RkCUyGiGTubp/uGH+fvzR+6Exxw/FDN6bVcyYZEXU+TIaICACaXXOIj8mIqDNjMkREAJreuLWsqhbF16oB2McGrURE5sZkiIgA3LTm0C09Q/m/vffqrkZ3jZPV4yIisjQmQ0QE4KaeoVvGDEkrT/MRGRF1UhZLhkpKSjB58mS4uLjAzc0NM2bMwLVrzS/sVlVVhbi4OPTs2RPdunXDhAkTUFhY2KDeunXrMGDAAGg0Gnh5eSEuLs5St0EkG8bB0YVl1bhWXSeVG3uK+nAbDiLqpCyWDE2ePBnHjx/Hrl278OWXX2LPnj2YNWtWs8fMnTsXX3zxBTZt2oRvvvkGFy9exOOPP25SZ8WKFViwYAHmz5+P48eP46uvvkJ0dLSlboNINlydneDRTQXgxqMx4Obd6tkzRESdk6MlTnry5ElkZGTg+++/x6BBgwAA77//PsaMGYPly5fDz8+vwTE6nQ5r1qzBhg0b8NBDDwEAPv74Y/Tt2xf79+/H0KFDcfXqVSxcuBBffPEFRo8eLR07YMAAS9wGkez08eiK4ms1yCu+hv69XAHc2MSVj8mIqLOySM9QdnY23NzcpEQIAKKiouDg4IADBw40ekxOTg5qa2sRFRUllYWGhiIwMBDZ2dkAgF27dsFgMODChQvo27cvevXqhSeffBLnz59vNp7q6mqUlZWZvIioIeMgauO6QgaDwNkrFSafERF1NhZJhrRaLby8vEzKHB0d0aNHD2i12iaPUalUcHNzMyn39vaWjsnLy4PBYMDSpUuRkpKCzZs3o6SkBL///e9RU1PTZDzJyclwdXWVXgEBAe27QaJO6tbp9ZfKqlBVa4CTUoFe7s62DI2IyGJalQzNnz8fCoWi2depU6csFSsMBgNqa2vx97//HdHR0Rg6dCg+/fRT/PLLL9i9e3eTxyUmJkKn00mv2/UkEcmVcRD1lz9dQtD87Ri+7GsAQO+eXeGo5ORTIuqcWjVm6KWXXsL06dObrRMcHAwfHx9cvnzZpLyurg4lJSXw8fFp9DgfHx/U1NSgtLTUpHeosLBQOsbX1xcAcM8990ife3p6wsPDAwUFBU3GpFaroVarm42biIDw3u7o2VWFKxWmPa1Rfb1tFBERkeW1Khny9PSEp6fnbetFRkaitLQUOTk5CA8PBwB8/fXXMBgMiIiIaPSY8PBwODk5ITMzExMmTAAA5ObmoqCgAJGRkQCA4cOHS+W9evUCUD+Fv7i4GL17927NrRBRI3p0VSE7cTTKq2qlMqWDAm5dVDaMiojIsizS7923b1888sgjmDlzJg4ePIjvvvsO8fHxeOqpp6SZZBcuXEBoaCgOHjwIAHB1dcWMGTOQkJCA3bt3IycnB7GxsYiMjMTQoUMBAHfddRcee+wxzJkzB/v27cOxY8cwbdo0hIaG4sEHH7TErRDJjsrRAT27qaUXEyEi6uwsNghg/fr1CA0NxejRozFmzBiMGDECq1atkj6vra1Fbm4uKisrpbK//e1v+OMf/4gJEybg/vvvh4+PD7Zs2WJy3k8++QQREREYO3YsRo0aBScnJ2RkZMDJidsEEMlNamoqgoKCoNFoEBERIf1y1ZRNmzYhNDQUGo0G/fv3x44dO6wUKRHZM4UQQtg6CGsrKyuDq6srdDodXFxcbB0OkSy193u4ceNGTJ06FWlpaYiIiEBKSgo2bdqE3NzcBrNZAWDfvn24//77kZycjD/+8Y/YsGED3n77bRw6dAj9+vWzSsxE1D6W+g4yGWKDRmQT7f0eRkREYPDgwVi5ciWA+tmmAQEBeOGFFzB//vwG9WNiYlBRUYEvv/xSKhs6dCjCwsKQlpZmlZiJqH0s9R3kXFki6nBqamqQk5Njskirg4MDoqKipEVab5WdnW1SHwCio6ObrE9E8mGR7TjsnbEzjCtRE9mO8fvXls7p4uJi6PV6eHubTvn39vZucq0zrVbbaP2mFoIF6levr66ult7rdDqT2InIutrTbjRHlslQeXk5AHAlaiI7UF5eDldXV1uH0ajk5GQsWbKkQTnbDiLbunLlilnbDVkmQ35+fjh//jy6d+8OhULRZL2ysjIEBATg/PnzHXZ8QEe/B8Zve5a6ByEEysvLG924+XY8PDygVCpRWFhoUn7zIq238vHxaVV9oH71+oSEBOl9aWkpevfujYKCArtN4G7VEf8OMmbr6Igx63Q6BAYGokePHmY9ryyTIQcHB2nRxpZwcXHpMH9RmtLR74Hx254l7qGtCYVKpUJ4eDgyMzMxfvx4APUDqDMzMxEfH9/oMZGRkcjMzMSLL74ole3atUta1LUxTa1e7+rq2uH+PDvi30HGbB0dMWYHB/MOeZZlMkREHV9CQgKmTZuGQYMGYciQIUhJSUFFRQViY2MBAFOnToW/vz+Sk5MBAHPmzMGoUaPw7rvvYuzYsUhPT8cPP/xgsv4ZEckTkyEi6pBiYmJQVFSExYsXQ6vVIiwsDBkZGdIg6YKCApPfHocNG4YNGzZg4cKFePXVV3HnnXdi27ZtLV5jiIg6LyZDzVCr1UhKSurQm7x29Htg/LZnz/cQHx/f5GOxrKysBmVPPPEEnnjiiTZfz55/Fk1hzNbBmK3DUjHLctFFIiIiIiMuukhERESyxmSIiIiIZI3JEBEREckakyEiIiKSNSZDzUhNTUVQUBA0Gg0iIiJw8OBBW4fUqOTkZAwePBjdu3eHl5cXxo8fj9zcXJM6DzzwABQKhclr9uzZNorY1GuvvdYgttDQUOnzqqoqxMXFoWfPnujWrRsmTJjQYCVhWwsKCmpwDwqFAnFxcQDs7+e/Z88ejBs3Dn5+flAoFNi2bZvJ50IILF68GL6+vnB2dkZUVBR++eUXkzolJSWYPHkyXFxc4ObmhhkzZuDatWtWvAvLaO33ftOmTQgNDYVGo0H//v2xY8cOK0V6Q2tiXr16NUaOHAl3d3e4u7sjKirKJm1bW9vX9PR0KBQKabFNa2ptzKWlpYiLi4Ovry/UajXuuusuq//9aG3MKSkpuPvuu+Hs7IyAgADMnTsXVVVVVor29m1TY7KysvC73/0OarUad9xxB9atW9f6CwtqVHp6ulCpVGLt2rXi+PHjYubMmcLNzU0UFhbaOrQGoqOjxccffyyOHTsmjhw5IsaMGSMCAwPFtWvXpDqjRo0SM2fOFJcuXZJeOp3OhlHfkJSUJO69916T2IqKiqTPZ8+eLQICAkRmZqb44YcfxNChQ8WwYcNsGHFDly9fNol/165dAoDYvXu3EML+fv47duwQCxYsEFu2bBEAxNatW00+X7ZsmXB1dRXbtm0TP/74o3j00UdFnz59xPXr16U6jzzyiBg4cKDYv3+/+Pbbb8Udd9whJk6caOU7Ma/Wfu+/++47oVQqxTvvvCNOnDghFi5cKJycnMTRo0ftNuZJkyaJ1NRUcfjwYXHy5Ekxffp04erqKn799Ve7jdkoPz9f+Pv7i5EjR4rHHnvMOsH+prUxV1dXi0GDBokxY8aIvXv3ivz8fJGVlSWOHDlitzGvX79eqNVqsX79epGfny927twpfH19xdy5c60W8+3aplvl5eWJLl26iISEBHHixAnx/vvvC6VSKTIyMlp1XSZDTRgyZIiIi4uT3uv1euHn5yeSk5NtGFXLXL58WQAQ33zzjVQ2atQoMWfOHNsF1YykpCQxcODARj8rLS0VTk5OYtOmTVLZyZMnBQCRnZ1tpQhbb86cOSIkJEQYDAYhhH3//G9tcAwGg/Dx8RF//etfpbLS0lKhVqvFp59+KoQQ4sSJEwKA+P7776U6//d//ycUCoW4cOGC1WI3t9Z+75988kkxduxYk7KIiAjx5z//2aJx3qy9bVVdXZ3o3r27+Mc//mGpEBtoS8x1dXVi2LBh4qOPPhLTpk2zejLU2pg//PBDERwcLGpqaqwVYgOtjTkuLk489NBDJmUJCQli+PDhFo2zKS1Jhl555RVx7733mpTFxMSI6OjoVl2Lj8kaUVNTg5ycHERFRUllDg4OiIqKQnZ2tg0jaxmdTgcADTayW79+PTw8PNCvXz8kJiaisrLSFuE16pdffoGfnx+Cg4MxefJkFBQUAABycnJQW1tr8mcRGhqKwMBAu/2zqKmpwb/+9S8888wzJhsB2/PP/2b5+fnQarUmP3NXV1dERERIP/Ps7Gy4ublh0KBBUp2oqCg4ODjgwIEDVo/ZHNryvc/OzjapDwDR0dFW+7tpjraqsrIStbW1Zt/4siltjfn111+Hl5cXZsyYYY0wTbQl5s8//xyRkZGIi4uDt7c3+vXrh6VLl0Kv19ttzMOGDUNOTo70KC0vLw87duzAmDFjrBJzW5jrO8gVqBtRXFwMvV4vLetv5O3tjVOnTtkoqpYxGAx48cUXMXz4cJNtBiZNmoTevXvDz88PP/30E+bNm4fc3Fxs2bLFhtHWi4iIwLp163D33Xfj0qVLWLJkCUaOHIljx45Bq9VCpVLBzc3N5Bhvb29otVrbBHwb27ZtQ2lpKaZPny6V2fPP/1bGn2tjf/+Nn2m1Wnh5eZl87ujoiB49etjtn8vttOV7r9Vqm/05WZo52qp58+bBz8+vwT8oltKWmPfu3Ys1a9bgyJEjVoiwobbEnJeXh6+//hqTJ0/Gjh07cPr0aTz//POora1FUlKSXcY8adIkFBcXY8SIERBCoK6uDrNnz8arr75q8XjbqqnvYFlZGa5fvw5nZ+cWnYfJUCcTFxeHY8eOYe/evSbls2bNkv6/f//+8PX1xejRo3HmzBmEhIRYO0wTf/jDH6T/HzBgACIiItC7d2989tlnLf6LbE/WrFmDP/zhD/Dz85PK7PnnT/K1bNkypKenIysrCxqNxtbhNKq8vBxTpkzB6tWr4eHhYetwWsxgMMDLywurVq2CUqlEeHg4Lly4gL/+9a9WSYbaIisrC0uXLsUHH3yAiIgInD59GnPmzMEbb7yBRYsW2To8i2Iy1AgPDw8olcoGM5YKCwvh4+Njo6huLz4+Hl9++SX27NmDXr16NVs3IiICAHD69Gm7+8fYzc0Nd911F06fPo3f//73qKmpQWlpqUnvkL3+WZw7dw5fffXVbXt87Pnnb/y5FhYWwtfXVyovLCxEWFiYVOfy5csmx9XV1aGkpMQu/1xaoi3fex8fH5u2E+1pq5YvX45ly5bhq6++woABAywZponWxnzmzBmcPXsW48aNk8oMBgOA+t7I3Nxci3+H2vJz9vX1hZOTE5RKpVTWt29faLVa1NTUQKVS2V3MixYtwpQpU/Dss88CqP/FraKiArNmzcKCBQtMNj62F019B11cXFr1y7T93ZkdUKlUCA8PR2ZmplRmMBiQmZmJyMhIG0bWOCEE4uPjsXXrVnz99dfo06fPbY8xdjff/I+dvbh27RrOnDkDX19fhIeHw8nJyeTPIjc3FwUFBXb5Z/Hxxx/Dy8sLY8eObbaePf/8+/TpAx8fH5OfeVlZGQ4cOCD9zCMjI1FaWoqcnBypztdffw2DwSAleh1NW773kZGRJvUBYNeuXVb7u9nWtuqdd97BG2+8gYyMDJNxX9bQ2phDQ0Nx9OhRHDlyRHo9+uijePDBB3HkyBEEBATYXcwAMHz4cJw+fVpK3ADg559/hq+vr8UTobbGXFlZ2SDhMSZzwk63MTXbd7BVw61lJD09XajVarFu3Tpx4sQJMWvWLOHm5ia0Wq2tQ2vgueeeE66uriIrK8tk6nZlZaUQQojTp0+L119/Xfzwww8iPz9f/Oc//xHBwcHi/vvvt3Hk9V566SWRlZUl8vPzxXfffSeioqKEh4eHuHz5shCifmp9YGCg+Prrr8UPP/wgIiMjRWRkpI2jbkiv14vAwEAxb948k3J7/PmXl5eLw4cPi8OHDwsAYsWKFeLw4cPi3LlzQoj6qfVubm7iP//5j/jpp5/EY4891ujU+vvuu08cOHBA7N27V9x5552dYmp9c9/7KVOmiPnz50v1v/vuO+Ho6CiWL18uTp48KZKSkmwytb41MS9btkyoVCqxefNmk/aivLzcbmO+lS1mk7U25oKCAtG9e3cRHx8vcnNzxZdffim8vLzEm2++abcxJyUlie7du4tPP/1U5OXlif/+978iJCREPPnkk1aL+XZt0/z588WUKVOk+sap9S+//LI4efKkSE1N5dR6c3v//fdFYGCgUKlUYsiQIWL//v22DqlRABp9ffzxx0KI+i/l/fffL3r06CHUarW44447xMsvv2w36wzFxMQIX19foVKphL+/v4iJiRGnT5+WPr9+/bp4/vnnhbu7u+jSpYv4n//5H3Hp0iUbRty4nTt3CgAiNzfXpNwef/67d+9u9O/MtGnThBD10+sXLVokvL29hVqtFqNHj25wX1euXBETJ04U3bp1Ey4uLiI2Ntaq/6BaSnPf+1GjRkk/I6PPPvtM3HXXXUKlUol7771XbN++3coRty7m3r17N/pnn5SUZLcx38oWyZAQrY953759IiIiQqjVahEcHCzeeustUVdXZ7cx19bWitdee02EhIQIjUYjAgICxPPPPy+uXr1qtXhv1zZNmzZNjBo1qsExYWFhQqVSieDgYOnfvtZQCGGnfV9EREREVsAxQ0RERCRrTIaIiIhI1pgMERERkawxGSIiIiJZYzJEREREssZkiIiIiGSNyRARERHJms2ToT179mDcuHHw8/ODQqHAtm3bbntMVlYWfve730GtVuOOO+7AunXrLB4nEdkPthtEZE42T4YqKiowcOBApKamtqh+fn4+xo4dK+1L8+KLL+LZZ5/Fzp07LRwpEdkLthtEZE52tQK1QqHA1q1bMX78+CbrzJs3D9u3b8exY8eksqeeegqlpaXIyMiwQpREZE/YbhBReznaOoDWys7ORlRUlElZdHQ0XnzxxSaPqa6uRnV1tfTeYDCgpKQEPXv2hEKhsFSoRNQMIQTKy8vh5+fXYKdsc2tLuwGw7SCyN5ZqNzpcMqTVauHt7W1S5u3tjbKyMly/fh3Ozs4NjklOTsaSJUusFSIRtcL58+fRq1cvi16jLe0GwLaDyF6Zu93ocMlQWyQmJiIhIUF6r9PpEBgYiPPnz8PFxcWGkRHJV1lZGQICAtC9e3dbh9Ikth1E9sVS7UaHS4Z8fHxQWFhoUlZYWAgXF5cmf7tTq9VQq9UNyl1cXNigEdmYNR43taXdANh2ENkrc7cbNp9N1lqRkZHIzMw0Kdu1axciIyNtFBER2Tu2G0TUHJsnQ9euXcORI0dw5MgRAPVTYI8cOYKCggIA9d3UU6dOlerPnj0beXl5eOWVV3Dq1Cl88MEH+OyzzzB37lxbhE9ENsB2g4jMStjY7t27BYAGr2nTpgkhhJg2bZoYNWpUg2PCwsKESqUSwcHB4uOPP27VNXU6nQAgdDqdeW6CiFqtPd9DW7Qb7Y2ZiNrPUt9Bu1pnyFrKysrg6uoKnU7H5/5ENtIRv4cdMWaizsRS30GbPyYjIiIisiUmQ0RERCRrTIaIiIhI1pgMERERkawxGSIiIiJZYzJEREREssZkiIiIiGSNyRARERHJGpMhIiIikjUmQ0RERCRrTIaIiIhI1pgMERERkawxGSIiIiJZYzJEREREssZkiIiIiGSNyRARERHJGpMhIiIikjUmQ0RERCRrTIaIiIhI1pgMERERkawxGSIiIiJZYzJEREREssZkiIiIiGTNLpKh1NRUBAUFQaPRICIiAgcPHmy2fkpKCu6++244OzsjICAAc+fORVVVlZWiJSJ7wbaDiMzB5snQxo0bkZCQgKSkJBw6dAgDBw5EdHQ0Ll++3Gj9DRs2YP78+UhKSsLJkyexZs0abNy4Ea+++qqVIyciW2LbQUTmYvNkaMWKFZg5cyZiY2Nxzz33IC0tDV26dMHatWsbrb9v3z4MHz4ckyZNQlBQEB5++GFMnDjxtr8RElHnwraDiMzFpslQTU0NcnJyEBUVJZU5ODggKioK2dnZjR4zbNgw5OTkSA1YXl4eduzYgTFjxjR5nerqapSVlZm8iKjjYttBRObkaMuLFxcXQ6/Xw9vb26Tc29sbp06davSYSZMmobi4GCNGjIAQAnV1dZg9e3azXd3JyclYsmSJWWMnItth20FE5mTzx2StlZWVhaVLl+KDDz7AoUOHsGXLFmzfvh1vvPFGk8ckJiZCp9NJr/Pnz1sxYiKyB2w7iKgpNu0Z8vDwgFKpRGFhoUl5YWEhfHx8Gj1m0aJFmDJlCp599lkAQP/+/VFRUYFZs2ZhwYIFcHBomN+p1Wqo1Wrz3wAR2QTbDiIyJ5v2DKlUKoSHhyMzM1MqMxgMyMzMRGRkZKPHVFZWNmi0lEolAEAIYblgichusO0gInOyac8QACQkJGDatGkYNGgQhgwZgpSUFFRUVCA2NhYAMHXqVPj7+yM5ORkAMG7cOKxYsQL33XcfIiIicPr0aSxatAjjxo2TGjYi6vzYdhCRudg8GYqJiUFRUREWL14MrVaLsLAwZGRkSAMjCwoKTH6bW7hwIRQKBRYuXIgLFy7A09MT48aNw1tvvWWrWyAiG2DbQUTmohAy7B8uKyuDq6srdDodXFxcbB0OkSx1xO9hR4yZqDOx1Heww80mIyIiIjInJkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLLGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLLGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLJmF8lQamoqgoKCoNFoEBERgYMHDzZbv7S0FHFxcfD19YVarcZdd92FHTt2WClaIrIXbDuIyBwcbR3Axo0bkZCQgLS0NERERCAlJQXR0dHIzc2Fl5dXg/o1NTX4/e9/Dy8vL2zevBn+/v44d+4c3NzcrB88EdkM2w4iMheFEELYMoCIiAgMHjwYK1euBAAYDAYEBATghRdewPz58xvUT0tLw1//+lecOnUKTk5ObbpmWVkZXF1dodPp4OLi0q74iaht2vs9ZNtBJD+W+g7a9DFZTU0NcnJyEBUVJZU5ODggKioK2dnZjR7z+eefIzIyEnFxcfD29ka/fv2wdOlS6PX6Jq9TXV2NsrIykxcRdVxsO4jInGyaDBUXF0Ov18Pb29uk3NvbG1qtttFj8vLysHnzZuj1euzYsQOLFi3Cu+++izfffLPJ6yQnJ8PV1VV6BQQEmPU+iMi62HYQkTnZxQDq1jAYDPDy8sKqVasQHh6OmJgYLFiwAGlpaU0ek5iYCJ1OJ73Onz9vxYiJyB6w7SCipth0ALWHhweUSiUKCwtNygsLC+Hj49PoMb6+vnBycoJSqZTK+vbtC61Wi5qaGqhUqgbHqNVqqNVq8wZPRDbDtoOIzMmmPUMqlQrh4eHIzMyUygwGAzIzMxEZGdnoMcOHD8fp06dhMBiksp9//hm+vr6NNmZE1Pmw7SAic7L5Y7KEhASsXr0a//jHP3Dy5Ek899xzqKioQGxsLABg6tSpSExMlOo/99xzKCkpwZw5c/Dzzz9j+/btWLp0KeLi4mx1C0RkA2w7iMhcbL7OUExMDIqKirB48WJotVqEhYUhIyNDGhhZUFAAB4cbOVtAQAB27tyJuXPnYsCAAfD398ecOXMwb948W90CEdkA2w4iMhebrzNkC1wrhMj2OuL3sCPGTNSZdMp1hoiIiIhsjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLLGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLLGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWbOLZCg1NRVBQUHQaDSIiIjAwYMHW3Rceno6FAoFxo8fb9kAicguse0gInOweTK0ceNGJCQkICkpCYcOHcLAgQMRHR2Ny5cvN3vc2bNn8Ze//AUjR460UqREZE/YdhCRudg8GVqxYgVmzpyJ2NhY3HPPPUhLS0OXLl2wdu3aJo/R6/WYPHkylixZguDgYCtGS0T2gm0HEZmLTZOhmpoa5OTkICoqSipzcHBAVFQUsrOzmzzu9ddfh5eXF2bMmNGi61RXV6OsrMzkRUQdF9sOIjInmyZDxcXF0Ov18Pb2Nin39vaGVqtt9Ji9e/dizZo1WL16dYuvk5ycDFdXV+kVEBDQrriJyLbYdhCROdn8MVlrlJeXY8qUKVi9ejU8PDxafFxiYiJ0Op30On/+vAWjJCJ7w7aDiJrjaMuLe3h4QKlUorCw0KS8sLAQPj4+DeqfOXMGZ8+exbhx46Qyg8EAAHB0dERubi5CQkIaHKdWq6FWq80cPRHZCtsOIjInm/YMqVQqhIeHIzMzUyozGAzIzMxEZGRkg/qhoaE4evQojhw5Ir0effRRPPjggzhy5Ai7sIlkgm0HEZmTTXuGACAhIQHTpk3DoEGDMGTIEKSkpKCiogKxsbEAgKlTp8Lf3x/JycnQaDTo16+fyfFubm4A0KCciDo3th1EZC42T4ZiYmJQVFSExYsXQ6vVIiwsDBkZGdLAyIKCAjg4dKihTURkBWw7iMhcFEIIYesgrK2srAyurq7Q6XRwcXGxdThEstQRv4cdMWaizsRS30H+2kRERESyxmSIiIiIZI3JEBEREckakyEiIiKSNSZDREREJGtMhoiIiEjWmAwRERGRrDEZIiIiIlljMkRERESyxmSIiIiIZI3JEBEREckakyEiIiKSNSZDREREJGtMhoiIiEjWmAwRERGRrDEZIiIiIlljMkRERESyxmSIiIiIZI3JEBEREckakyEiIiKSNSZDREREJGtMhoiIiEjW7CIZSk1NRVBQEDQaDSIiInDw4MEm665evRojR46Eu7s73N3dERUV1Wx9Iuq82HYQkTnYPBnauHEjEhISkJSUhEOHDmHgwIGIjo7G5cuXG62flZWFiRMnYvfu3cjOzkZAQAAefvhhXLhwwcqRE5Etse0gInNRCCGELQOIiIjA4MGDsXLlSgCAwWBAQEAAXnjhBcyfP/+2x+v1eri7u2PlypWYOnVqi65ZVlYGV1dX6HQ6uLi4tCt+Imqb9n4P2XYQyY+lvoM27RmqqalBTk4OoqKipDIHBwdERUUhOzu7ReeorKxEbW0tevToYakwicjOsO0gInNytOXFi4uLodfr4e3tbVLu7e2NU6dOtegc8+bNg5+fn0mjeKvq6mpUV1dL78vKytoWMBHZBbYdRGRONh8z1B7Lli1Deno6tm7dCo1G02S95ORkuLq6Sq+AgAArRklE9oZtBxHdzKbJkIeHB5RKJQoLC03KCwsL4ePj0+yxy5cvx7Jly/Df//4XAwYMaLZuYmIidDqd9Dp//ny7Yyci22HbQUTmZNNkSKVSITw8HJmZmVKZwWBAZmYmIiMjmzzunXfewRtvvIGMjAwMGjTottdRq9VwcXExeRFRx8W2g4jMyaZjhgAgISEB06ZNw6BBgzBkyBCkpKSgoqICsbGxAICpU6fC398fycnJAIC3334bixcvxoYNGxAUFAStVgsA6NatG7p162az+yAi62LbQUTmYvNkKCYmBkVFRVi8eDG0Wi3CwsKQkZEhDYwsKCiAg8ONDqwPP/wQNTU1+NOf/mRynqSkJLz22mvWDJ2IbIhtBxGZi83XGbIFrhVCZHsd8XvYEWMm6kw65TpDRERERLbGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLLGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGsMRkiIiIiWWMyRERERLLGZIiIiIhkjckQERERyRqTISIiIpI1JkNEREQka0yGiIiISNaYDBEREZGs2UUylJqaiqCgIGg0GkRERODgwYPN1t+0aRNCQ0Oh0WjQv39/7Nixw0qREpE9YdtBROZg82Ro48aNSEhIQFJSEg4dOoSBAwciOjoaly9fbrT+vn37MHHiRMyYMQOHDx/G+PHjMX78eBw7dszKkRORLbHtICJzUQghhC0DiIiIwODBg7Fy5UoAgMFgQEBAAF544QXMnz+/Qf2YmBhUVFTgyy+/lMqGDh2KsLAwpKWlteiaZWVlcHV1hU6ng4uLi3luhIhapb3fQ7YdRPJjqe+go9nO1AY1NTXIyclBYmKiVObg4ICoqChkZ2c3ekx2djYSEhJMyqKjo7Ft27Ymr1NdXY3q6mrpvU6nA1D/QyUi2zB+/9ry+xjbDiJ5ak+70RybJkPFxcXQ6/Xw9vY2Kff29sapU6caPUar1TZaX6vVNnmd5ORkLFmypEF5QEBAG6ImInO6cuUKXF1dW3UM2w4ieWtLu9EcmyZD1pKYmGjyG2FpaSl69+6NgoICs/4wLamsrAwBAQE4f/58h+meZ8zW0RFjBup7WQIDA9GjRw9bh9Ikth22wZitoyPGbKl2w6bJkIeHB5RKJQoLC03KCwsL4ePj0+gxPj4+raoPAGq1Gmq1ukG5q6trh/kLYOTi4sKYrYAxW4+DQ+vncbDtaL2O+PeDMVtHR4y5Le1Gs+cz69laSaVSITw8HJmZmVKZwWBAZmYmIiMjGz0mMjLSpD4A7Nq1q8n6RNT5sO0gInOy+WOyhIQETJs2DYMGDcKQIUOQkpKCiooKxMbGAgCmTp0Kf39/JCcnAwDmzJmDUaNG4d1338XYsWORnp6OH374AatWrbLlbRCRlbHtICJzsXkyFBMTg6KiIixevBharRZhYWHIyMiQBjoWFBSYdIcNGzYMGzZswMKFC/Hqq6/izjvvxLZt29CvX78WX1OtViMpKanR7m97xZitgzFbT3vjZtvRMozZOhizdVgqZpuvM0RERERkSzZfgZqIiIjIlpgMERERkawxGSIiIiJZYzJEREREstZpk6HU1FQEBQVBo9EgIiICBw8ebLb+pk2bEBoaCo1Gg/79+2PHjh1WivSG1sS8evVqjBw5Eu7u7nB3d0dUVNRt79ESWvtzNkpPT4dCocD48eMtG2AjWhtzaWkp4uLi4OvrC7Vajbvuusvqfz9aG3NKSgruvvtuODs7IyAgAHPnzkVVVZWVogX27NmDcePGwc/PDwqFotn9v4yysrLwu9/9Dmq1GnfccQfWrVtn8Tgbw7bDOth2WAfbjhYSnVB6erpQqVRi7dq14vjx42LmzJnCzc1NFBYWNlr/u+++E0qlUrzzzjvixIkTYuHChcLJyUkcPXrUbmOeNGmSSE1NFYcPHxYnT54U06dPF66uruLXX3+125iN8vPzhb+/vxg5cqR47LHHrBPsb1obc3V1tRg0aJAYM2aM2Lt3r8jPzxdZWVniyJEjdhvz+vXrhVqtFuvXrxf5+fli586dwtfXV8ydO9dqMe/YsUMsWLBAbNmyRQAQW7dubbZ+Xl6e6NKli0hISBAnTpwQ77//vlAqlSIjI8M6Af+GbYd9xmzEtsOyMcu57eiUydCQIUNEXFyc9F6v1ws/Pz+RnJzcaP0nn3xSjB071qQsIiJC/PnPf7ZonDdrbcy3qqurE927dxf/+Mc/LBViA22Jua6uTgwbNkx89NFHYtq0aVZv0Fob84cffiiCg4NFTU2NtUJsoLUxx8XFiYceesikLCEhQQwfPtyicTalJQ3aK6+8Iu69916TspiYGBEdHW3ByBpi22EdbDusg21Hy3W6x2Q1NTXIyclBVFSUVObg4ICoqChkZ2c3ekx2drZJfQCIjo5usr65tSXmW1VWVqK2ttZqm162NebXX38dXl5emDFjhjXCNNGWmD///HNERkYiLi4O3t7e6NevH5YuXQq9Xm+3MQ8bNgw5OTlSd3heXh527NiBMWPGWCXmtrD1dxBg28G2o2lsOzp/22HzFajNrbi4GHq9XlqF1sjb2xunTp1q9BitVttofa1Wa7E4b9aWmG81b948+Pn5NfhLYSltiXnv3r1Ys2YNjhw5YoUIG2pLzHl5efj6668xefJk7NixA6dPn8bzzz+P2tpaJCUl2WXMkyZNQnFxMUaMGAEhBOrq6jB79my8+uqrFo+3rZr6DpaVleH69etwdna2eAxsO9h2NIVtR+dvOzpdz5AcLVu2DOnp6di6dSs0Go2tw2lUeXk5pkyZgtWrV8PDw8PW4bSYwWCAl5cXVq1ahfDwcMTExGDBggVIS0uzdWhNysrKwtKlS/HBBx/g0KFD2LJlC7Zv34433njD1qGRnWHbYTlsOzqWTtcz5OHhAaVSicLCQpPywsJC+Pj4NHqMj49Pq+qbW1tiNlq+fDmWLVuGr776CgMGDLBkmCZaG/OZM2dw9uxZjBs3TiozGAwAAEdHR+Tm5iIkJMSuYgYAX19fODk5QalUSmV9+/aFVqtFTU0NVCqV3cW8aNEiTJkyBc8++ywAoH///qioqMCsWbOwYMECk/267EVT30EXFxer9AoBbDushW0H2w5zMlfbYX931k4qlQrh4eHIzMyUygwGAzIzMxEZGdnoMZGRkSb1AWDXrl1N1je3tsQMAO+88w7eeOMNZGRkYNCgQdYIVdLamENDQ3H06FEcOXJEej366KN48MEHceTIEQQEBNhdzAAwfPhwnD59Wmp8AeDnn3+Gr6+vxRuztsZcWVnZoNEyNsjCTrcitPV3EGDbYS1sO9h2mJPZvoOtGm7dQaSnpwu1Wi3WrVsnTpw4IWbNmiXc3NyEVqsVQggxZcoUMX/+fKn+d999JxwdHcXy5cvFyZMnRVJSkk2mx7Ym5mXLlgmVSiU2b94sLl26JL3Ky8vtNuZb2WJGSGtjLigoEN27dxfx8fEiNzdXfPnll8LLy0u8+eabdhtzUlKS6N69u/j0009FXl6e+O9//ytCQkLEk08+abWYy8vLxeHDh8Xhw4cFALFixQpx+PBhce7cOSGEEPPnzxdTpkyR6hunx7788svi5MmTIjU11WZT69l22F/Mt2LbYZmY5dx2dMpkSAgh3n//fREYGChUKpUYMmSI2L9/v/TZqFGjxLRp00zqf/bZZ+Kuu+4SKpVK3HvvvWL79u1Wjrh1Mffu3VsAaPBKSkqy25hvZYsGTYjWx7xv3z4REREh1Gq1CA4OFm+99Zaoq6uz25hra2vFa6+9JkJCQoRGoxEBAQHi+eefF1evXrVavLt3727076cxzmnTpolRo0Y1OCYsLEyoVCoRHBwsPv74Y6vFezO2HfYX863YdlgmZjm3HQoh7LTvi4iIiMgKOt2YISIiIqLWYDJEREREssZkiIiIiGSNyRARERHJGpMhIiIikjUmQ0RERCRrTIaIiIhI1pgMERERkawxGSIiIiJZYzJEREREssZkiIiIiGSNyRARERHJ2v8HebDUhQBG3RIAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import csv\n", + "from math import pi\n", + "import matplotlib.pyplot as plt\n", + "\n", + "angle = []\n", + "\n", + "with open('../log/11-30_05%3A35%3A26new.csv', newline='') as csvfile:\n", + " reader = csv.reader(csvfile, delimiter=',')\n", + " for row in reader:\n", + " if row[1] != \"Angle\":\n", + " angle.append(float(row[1]) * (pi/360))\n", + " \n", + " fig, axs = plt.subplots(2, 2)\n", + " fig.suptitle(\"Run 1\")\n", + " \n", + " axs[0,0].plot(angle)\n", + " axs[0,0].set_title('Angle [rad]')\n", + " # axs[0,1].plot(self.v_ang)\n", + " # axs[0,1].set_title('Angular velocity [rad/s]')\n", + " # axs[1,0].plot(self.a_ang)\n", + " # axs[1,0].set_title('Angular acceleration [rad/s^2]')\n", + " # axs[1,1].plot(self.a_ang)\n", + " # axs[1,1].set_title('Angular acceleration [rad/s^2]')\n", + " " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "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.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/res/bricklink/modcon_balancer.io b/res/bricklink/modcon_balancer.io index 76239fc..66133ad 100644 Binary files a/res/bricklink/modcon_balancer.io and b/res/bricklink/modcon_balancer.io differ diff --git a/res/drawio/system_diagram.drawio b/res/drawio/system_diagram.drawio new file mode 100644 index 0000000..c8c1bad --- /dev/null +++ b/res/drawio/system_diagram.drawio @@ -0,0 +1,100 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/res/drawio/system_diagram.png b/res/drawio/system_diagram.png new file mode 100644 index 0000000..0e35a67 Binary files /dev/null and b/res/drawio/system_diagram.png differ diff --git a/res/screenshot2.png b/res/screenshot2.png new file mode 100644 index 0000000..f648413 Binary files /dev/null and b/res/screenshot2.png differ diff --git a/src/ev3.py b/src/ev3.py deleted file mode 100644 index f6506c2..0000000 --- a/src/ev3.py +++ /dev/null @@ -1,26 +0,0 @@ -#!/usr/bin/env python3 -import os -import time - -from ev3dev2.motor import LargeMotor, OUTPUT_B, OUTPUT_C, SpeedPercent, MoveTank -from ev3dev2.sensor import INPUT_2 -from ev3dev2.sensor.lego import GyroSensor -from ev3dev2.led import Leds - -os.system('setfont Lat15-TerminusBold14') - -# tank_drive = MoveTank(OUTPUT_B,OUTPUT_C) - -gyro = GyroSensor(INPUT_2) - -dt = 0 -t = time.time() - -while (True): - # Timekeepers - dt = time.time() - t - t = time.time() - - # Gyro - gyro_value = gyro.angle - print("Gyro: " + str(gyro_value)) diff --git a/src/ev3/CMakeLists.txt b/src/ev3_c/CMakeLists.txt similarity index 100% rename from src/ev3/CMakeLists.txt rename to src/ev3_c/CMakeLists.txt diff --git a/src/ev3/cc.ps1 b/src/ev3_c/cc.ps1 similarity index 100% rename from src/ev3/cc.ps1 rename to src/ev3_c/cc.ps1 diff --git a/src/ev3/cli.ps1 b/src/ev3_c/cli.ps1 similarity index 100% rename from src/ev3/cli.ps1 rename to src/ev3_c/cli.ps1 diff --git a/src/ev3/deploy.ps1 b/src/ev3_c/deploy.ps1 similarity index 100% rename from src/ev3/deploy.ps1 rename to src/ev3_c/deploy.ps1 diff --git a/doc/ethics.tex b/src/ev3_c/src/delay.c similarity index 100% rename from doc/ethics.tex rename to src/ev3_c/src/delay.c diff --git a/src/ev3/src/main.c b/src/ev3_c/src/main.c similarity index 100% rename from src/ev3/src/main.c rename to src/ev3_c/src/main.c diff --git a/src/ev3_py/ev3.py b/src/ev3_py/ev3.py new file mode 100644 index 0000000..35757bc --- /dev/null +++ b/src/ev3_py/ev3.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python3 +import os +import time +import sys + +from ev3dev2.motor import LargeMotor, OUTPUT_B, OUTPUT_C, SpeedPercent, MoveTank +from ev3dev2.sensor import INPUT_2 +from ev3dev2.sensor.lego import GyroSensor +from ev3dev2.button import Button + +import time +import logging +from datetime import datetime + +now = datetime.now() # current date and time +LOGNAME = os.path.join(os.getcwd(), "logs/", (now.strftime("%m-%d_%H:%M:%S") + ".csv")) +# print(LOGNAME) +logging.basicConfig(filename=LOGNAME, + filemode='a', + format='%(message)s', + level=logging.DEBUG) + +os.system('setfont Lat15-TerminusBold14') + +# Motors +tank_drive = MoveTank(OUTPUT_B,OUTPUT_C) + +def exit(): + print('Adios') + tank_drive.stop() + sys.exit(130) + +# Button +btn = Button() +btn.on_backspace = exit + +# Gyro +gyro = GyroSensor(INPUT_2) +gyro.mode = GyroSensor.MODE_GYRO_RATE + +print("Put system in stable position!") +for i in range(1, 0, -1): + print(i) + time.sleep(1) + +gyro.calibrate() + +# Keep time +dt = 0 +t = time.time() +error = 0 + +def main(): + global dt, t, error + + Kp = 3 + Ki = 0 + Kd = 0.15 + + power = 0 + while (True): + # Timekeepers + dt = time.time() - t + t = time.time() + + # Gyro + prev_err = error + gyro_value = gyro.rate + error += (gyro_value * dt) + + power = (error * Kp) + (((error - prev_err) / (dt/1000)) * Kd) + # logging.debug("%d;%d;%d", error, power) + # print("Err: ", error, "Pwr: ", power) + + if (power < -100): + tank_drive.on(SpeedPercent(-100), SpeedPercent(-100)) + elif (power > 100): + tank_drive.on(SpeedPercent(100), SpeedPercent(100)) + else: + tank_drive.on(SpeedPercent(power), SpeedPercent(power)) + + # Has likely fallen over + if (error > 40 or error < -40): + exit() + +if __name__ == '__main__': + try: + main() + except KeyboardInterrupt: + exit() \ No newline at end of file diff --git a/src/ev3_py/getGyro.py b/src/ev3_py/getGyro.py new file mode 100644 index 0000000..83ba163 --- /dev/null +++ b/src/ev3_py/getGyro.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python3 +import os +import time + +from ev3dev2.motor import LargeMotor, OUTPUT_B, OUTPUT_C, SpeedPercent, MoveTank +from ev3dev2.sensor import INPUT_2 +from ev3dev2.sensor.lego import GyroSensor +from ev3dev2.led import Leds + +os.system('setfont Lat15-TerminusBold14') + +# tank_drive = MoveTank(OUTPUT_B,OUTPUT_C) + +gyro = GyroSensor(INPUT_2) +gyro.mode = GyroSensor.MODE_GYRO_RATE + +dt = 0 +t = time.time() + +error = 0.0 +total = 0 +for i in range(10): + total += gyro.rate + time.sleep(0.01) +average = total / 10.0 +offset = average - 0.1 + +GYRO_OFFSET = 0.0005 + +while (True): + # Timekeepers + dt = time.time() - t + t = time.time() + + gyro_raw = gyro.rate + offset = GYRO_OFFSET * gyro_raw + (1 - GYRO_OFFSET) * offset + speed = (gyro_raw - offset) + error += speed * (dt/1000) + + # # Gyro + # gyro_value = gyro.rate + # error += (gyro_value * dt) + + print("Gyro: " + str(error)) \ No newline at end of file diff --git a/src/ev3_py/helpers.py b/src/ev3_py/helpers.py new file mode 100644 index 0000000..ac82e6d --- /dev/null +++ b/src/ev3_py/helpers.py @@ -0,0 +1,4 @@ +def exit(): + print('Adios') + tank_drive.stop() + sys.exit(130) diff --git a/src/htway.py b/src/ev3_py/htway.py similarity index 95% rename from src/htway.py rename to src/ev3_py/htway.py index ee8c99e..ac4d8d2 100644 --- a/src/htway.py +++ b/src/ev3_py/htway.py @@ -3,6 +3,8 @@ # Python port of the HiTechnic HTWay for ev3dev # Copyright (c) 2014-2015 G33kDude, David Lechner # HiTechnic HTWay is Copyright (c) 2010 HiTechnic + + import itertools import os import struct @@ -22,14 +24,14 @@ WHEEL_RATIO_RCX = 1.4 # 81.6mm # the wheel size. KGYROANGLE = 7.5 KGYROSPEED = 1.15 -KPOS = 0.07 -KSPEED = 0.1 +KPOS = 0 +KSPEED = 0 # This constant aids in drive control. When the robot starts moving # because of user control, this constant helps get the robot leaning in # the right direction. Similarly, it helps bring robot to a stop when # stopping. -KDRIVE = -0.02 +KDRIVE = 0 # Power differential used for steering based on difference of target # steering and actual motor difference. @@ -84,6 +86,7 @@ class Gyro: self.value0.seek(0) return int(self.value0.read()) + class EV3Gyro(Gyro): DRIVER_NAME = 'lego-ev3-gyro' @@ -111,23 +114,27 @@ class HTGyro(Gyro): self.angle += speed * interval return speed, self.angle + class EV3Motor: def __init__(self, which=0): - devices = list(pyudev.Context().list_devices(subsystem='tacho-motor').match_attribute('driver_name', 'lego-ev3-l-motor')) + devices = list(pyudev.Context().list_devices(subsystem='tacho-motor') \ + .match_attribute('driver_name', 'lego-ev3-l-motor')) if which >= len(devices): raise Exception("Motor not found") self.device = devices[which] - self.pos = open(os.path.join(self.device.sys_path, 'position'), 'r+',0) - self.duty_cycle_sp = open(os.path.join(self.device.sys_path,'duty_cycle_sp'), 'w', 0) + self.pos = open(os.path.join(self.device.sys_path, 'position'), 'r+', + 0) + self.duty_cycle_sp = open(os.path.join(self.device.sys_path, + 'duty_cycle_sp'), 'w', 0) + self.reset() def reset(self): self.set_pos(0) self.set_duty_cycle_sp(0) self.send_command("run-direct") - def get_pos(self): self.pos.seek(0) return int(self.pos.read()) @@ -139,7 +146,8 @@ class EV3Motor: return self.duty_cycle_sp.write(str(int(new_duty_cycle_sp))) def send_command(self, new_mode): - with open(os.path.join(self.device.sys_path, 'command'),'w') as command: + with open(os.path.join(self.device.sys_path, 'command'), + 'w') as command: command.write(str(new_mode)) class EV3Motors: @@ -185,23 +193,18 @@ def main(): gyro = Gyro.get_gyro() gyro.calibrate() - - print("balancing in ...") + print "balancing in ..." for i in range(5, 0, -1): - print(i) + print i os.system("beep -f 440 -l 100") time.sleep(1) - print(0) - + print 0 os.system("beep -f 440 -l 1000") time.sleep(1) - motors = EV3Motors() start_time = time.time() last_okay_time = start_time - avg_interval = 0.0055 - for loop_count in itertools.count(): gyro_speed, gyro_angle = gyro.get_data(avg_interval) motors_speed, motors_pos = motors.get_data(avg_interval) @@ -225,6 +228,7 @@ def main(): time.sleep(WAIT_TIME) + now_time = time.time() avg_interval = (now_time - start_time) / (loop_count+1) @@ -237,4 +241,4 @@ def main(): motors.right.send_command('reset') if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/src/ev3_py/htway_og.py b/src/ev3_py/htway_og.py new file mode 100644 index 0000000..de69779 --- /dev/null +++ b/src/ev3_py/htway_og.py @@ -0,0 +1,244 @@ +#!/usr/bin/env python + +# Python port of the HiTechnic HTWay for ev3dev +# Copyright (c) 2014-2015 G33kDude, David Lechner +# HiTechnic HTWay is Copyright (c) 2010 HiTechnic + + +import itertools +import os +import struct +import time +import pyudev + +# loop interval in seconds +WAIT_TIME = 0.008 + +# ratio of wheel size compared to NXT 1.0 +WHEEL_RATIO_NXT1 = 1.0 # 56mm +WHEEL_RATIO_NXT = 0.8 # 43.2mm (same as EV3) +WHEEL_RATIO_RCX = 1.4 # 81.6mm + +# These are the main four balance constants, only the gyro constants +# are relative to the wheel size. KPOS and KSPEED are self-relative to +# the wheel size. +KGYROANGLE = 7.5 +KGYROSPEED = 1.15 +KPOS = 0.07 +KSPEED = 0.1 + +# This constant aids in drive control. When the robot starts moving +# because of user control, this constant helps get the robot leaning in +# the right direction. Similarly, it helps bring robot to a stop when +# stopping. +KDRIVE = -0.02 + +# Power differential used for steering based on difference of target +# steering and actual motor difference. +KSTEER = 0.25 + +# If robot power is saturated (over +/- 100) for over this time limit +# then robot must have fallen. In seconds. +TIME_FALL_LIMIT = 0.5 + +# Gyro offset control +# The HiTechnic gyro sensor will drift with time. This constant is +# used in a simple long term averaging to adjust for this drift. +# Every time through the loop, the current gyro sensor value is +# averaged into the gyro offset weighted according to this constant. +EMAOFFSET = 0.0005 + +class Gyro: + @staticmethod + def get_gyro(): + devices = list(pyudev.Context().list_devices(subsystem='lego-sensor') \ + .match_property("LEGO_DRIVER_NAME", EV3Gyro.DRIVER_NAME) \ + .match_property("LEGO_DRIVER_NAME", HTGyro.DRIVER_NAME)) + + if not devices: + raise Exception('Gyro not found') + + if devices[0].attributes.get('driver_name') == EV3Gyro.DRIVER_NAME: + return EV3Gyro(devices[0]) + elif devices[0].attributes.get('driver_name') == HTGyro.DRIVER_NAME: + return HTGyro(devices[0]) + + def __init__(self, device): + self.device = device + self.value0 = open(os.path.join(self.device.sys_path, 'value0'), 'r', 0) + self.offset = 0.0 + self.angle = 0.0 + + def calibrate(self): + self.angle = 0.0 + total = 0 + for i in range(10): + total += self.get_rate() + time.sleep(0.01) + average = total / 10.0 + self.offset = average - 0.1 + + def set_mode(self, value): + with open(os.path.join(self.device.sys_path, 'mode'), 'w') as f: + f.write(str(value)) + + def get_rate(self): + self.value0.seek(0) + return int(self.value0.read()) + + +class EV3Gyro(Gyro): + DRIVER_NAME = 'lego-ev3-gyro' + + def __init__(self, device): + Gyro.__init__(self, device) + self.set_mode("GYRO-RATE") + + def get_data(self, interval): + gyro_raw = self.get_rate() + self.offset = EMAOFFSET * gyro_raw + (1 - EMAOFFSET) * self.offset + speed = (gyro_raw - self.offset) + self.angle += speed * interval + return speed, self.angle + +class HTGyro(Gyro): + DRIVER_NAME = 'ht-nxt-gyro' + + def __init__(self, device): + Gyro.__init__(self, device) + + def get_data(self, interval): + gyro_raw = self.get_rate() + self.offset = EMAOFFSET * gyro_raw + (1 - EMAOFFSET) * self.offset + speed = (gyro_raw - self.offset) * 400.0 / 1953.0 + self.angle += speed * interval + return speed, self.angle + + +class EV3Motor: + def __init__(self, which=0): + devices = list(pyudev.Context().list_devices(subsystem='tacho-motor') \ + .match_attribute('driver_name', 'lego-ev3-l-motor')) + + if which >= len(devices): + raise Exception("Motor not found") + + self.device = devices[which] + self.pos = open(os.path.join(self.device.sys_path, 'position'), 'r+', + 0) + self.duty_cycle_sp = open(os.path.join(self.device.sys_path, + 'duty_cycle_sp'), 'w', 0) + + self.reset() + + def reset(self): + self.set_pos(0) + self.set_duty_cycle_sp(0) + self.send_command("run-direct") + def get_pos(self): + self.pos.seek(0) + return int(self.pos.read()) + + def set_pos(self, new_pos): + return self.pos.write(str(int(new_pos))) + + def set_duty_cycle_sp(self, new_duty_cycle_sp): + return self.duty_cycle_sp.write(str(int(new_duty_cycle_sp))) + + def send_command(self, new_mode): + with open(os.path.join(self.device.sys_path, 'command'), + 'w') as command: + command.write(str(new_mode)) + +class EV3Motors: + def __init__(self, left=0, right=1): + self.left = EV3Motor(left) + self.right = EV3Motor(right) + self.pos = 0.0 + self.speed = 0.0 + self.diff = 0 + self.target_diff = 0 + self.drive = 0 + self.steer = 0 + self.prev_sum = 0 + self.prev_deltas = [0,0,0] + + def get_data(self, interval): + left_pos = self.left.get_pos() + right_pos = self.right.get_pos() + + pos_sum = right_pos + left_pos + self.diff = left_pos - right_pos + + delta = pos_sum - self.prev_sum + self.pos += delta + + self.speed = (delta + sum(self.prev_deltas)) / (4 * interval) + + self.prev_sum = pos_sum + self.prev_deltas = [delta] + self.prev_deltas[0:2] + + return self.speed, self.pos + + def steer_control(self, power, steering, interval): + self.target_diff += steering * interval + power_steer = KSTEER * (self.target_diff - self.diff) + power_left = max(-100, min(power + power_steer, 100)) + power_right = max(-100, min(power - power_steer, 100)) + return power_left, power_right + + +def main(): + wheel_ratio = WHEEL_RATIO_NXT + + gyro = Gyro.get_gyro() + gyro.calibrate() + print "balancing in ..." + for i in range(5, 0, -1): + print i + os.system("beep -f 440 -l 100") + time.sleep(1) + print 0 + os.system("beep -f 440 -l 1000") + time.sleep(1) + motors = EV3Motors() + start_time = time.time() + last_okay_time = start_time + avg_interval = 0.0055 + for loop_count in itertools.count(): + gyro_speed, gyro_angle = gyro.get_data(avg_interval) + motors_speed, motors_pos = motors.get_data(avg_interval) + #print gyro_speed, gyro_angle, motors_speed, motors_pos + motors_pos -= motors.drive * avg_interval + motors.pos = motors_pos + + power = (KGYROSPEED * gyro_speed + + KGYROANGLE * gyro_angle) \ + / wheel_ratio \ + + KPOS * motors_pos \ + + KDRIVE * motors.drive \ + + KSPEED * motors_speed + + left_power, right_power = motors.steer_control(power, 0, avg_interval) + + #print left_power, right_power + + motors.left.set_duty_cycle_sp(left_power) + motors.right.set_duty_cycle_sp(right_power) + + time.sleep(WAIT_TIME) + + + now_time = time.time() + avg_interval = (now_time - start_time) / (loop_count+1) + + if abs(power) < 100: + last_okay_time = now_time + elif now_time - last_okay_time > TIME_FALL_LIMIT: + break + + motors.left.send_command('reset') + motors.right.send_command('reset') + +if __name__ == "__main__": + main() diff --git a/src/ev3_py/test.py b/src/ev3_py/test.py new file mode 100644 index 0000000..141c09a --- /dev/null +++ b/src/ev3_py/test.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 + +# https://github.com/ev3dev/ev3dev-lang-python +# https://github.com/ev3dev/ev3dev-lang-python-demo#balanc3r + +from ev3dev2.motor import LargeMotor, OUTPUT_B, OUTPUT_C +from ev3dev2.sensor import INPUT_2 +from ev3dev2.sensor.lego import GyroSensor +from ev3dev2.sound import Sound + +import time +# import logging +import os +import csv + +from datetime import datetime + +now = datetime.now() # current date and time +# LOGNAME = os.path.join(os.getcwd(), "logs/", (now.strftime("%m-%d_%H:%M:%S") + "final.csv")) +# # print(LOGNAME) +# logging.basicConfig(filename=LOGNAME, +# filemode='a', +# format='%(message)s', +# level=logging.DEBUG) + +log_t = [] +log_ang = [] +log_vang = [] +log_x = [] +log_v = [] + +# Motors +l_motor = LargeMotor(OUTPUT_B) +r_motor = LargeMotor(OUTPUT_C) + +# Sensors +gyro = GyroSensor(INPUT_2) + +# Sound +sound = Sound() + +# tires are 56mm diameter +def balance(target_angle=0): + t = time.time() + + # PID + gyro_k_p = 7.5 + # float k_i = 0 # No integral in this system + gyro_k_d = 1.15 + motor_k_p = 0.07 + motor_k_d = 0.1 + + # Calibrate gyro in current position + gyro.calibrate() + gyro.mode = GyroSensor.MODE_GYRO_G_A + angle = gyro.angle + + # Calibrate motor + prev_sum_motor_pos = 0 + l_motor.reset() + r_motor.reset() + pos = 0 + + sound.speak('3, 2, 1') + + # Stop if the robot has fallen over + while abs(angle) < 40: + # Keep time + dt = time.time() - t + t = time.time() + + sum_motor_pos = l_motor.position + r_motor.position + + delta_motor_pos = sum_motor_pos - prev_sum_motor_pos + pos += delta_motor_pos + speed = (delta_motor_pos / dt) + + prev_sum_motor_pos = sum_motor_pos + + angle, rate = gyro.angle_and_rate + + log_t.append(t) + log_ang.append(angle) + log_vang.append(rate) + log_x.append(pos) + log_v.append(speed) + + error = angle - target_angle + pd = gyro_k_p * error + gyro_k_d * rate \ + + motor_k_p * pos + motor_k_d * speed + + # convert -100 ~ 100 to -1050 ~ 1050 + speed = (((pd - (-100)) * (1049 - (-1049))) / (100 - (-100))) + (-1049) + l_motor.run_forever(speed_sp=speed) + r_motor.run_forever(speed_sp=speed) + + # logging.debug("%d;%d;%d", angle, rate, pd) + + # if abs(pd) > 1050: + # print("cap") + + # speed = min(max(pd, -1050), 1050) + + l_motor.stop(stop_action="hold") + r_motor.stop(stop_action="hold") + + +# Combine the arrays +try: + balance() +except: + l_motor.stop(stop_action="hold") + r_motor.stop(stop_action="hold") +finally: + data = zip(log_t, log_ang, log_vang, log_x, log_v) + + # Define the filename + now = datetime.now() # current date and time + LOGNAME = os.path.join(now.strftime("%m-%d_%H:%M:%S") + "new.csv") + + # Write data to CSV file + with open(LOGNAME, 'w', newline='') as csvfile: + csvwriter = csv.writer(csvfile) + csvwriter.writerow(['Time', 'Angle', 'Angular Velocity', 'X pos', 'Speed']) # Write header row + csvwriter.writerows(data) diff --git a/src/gyro.py b/src/gyro.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/sim/pendulum.py b/src/sim/pendulum.py index 5eaee43..cc255e5 100644 --- a/src/sim/pendulum.py +++ b/src/sim/pendulum.py @@ -10,6 +10,9 @@ C_MTPRATIO = 100 # Pixels per meter C_P_ANG_START = 1 / 1000 * math.pi C_FALL_ANG = 52.5 / 100 * math.pi +C_ANG_RES_LIMIT = ((2 * math.pi) / 360) # rad +C_V_CART_LIMIT = 0.403 # m/s +C_A_CART_LIMIT = C_V_CART_LIMIT / 0.1 # m/s^2 class Pendulum: def __init__(self, theta, length, dx, mass, color): @@ -35,6 +38,9 @@ class Pendulum: self.theta = [theta] # Angle in radians self.a_ang = [0] # Angular acceleration self.v_ang = [0] # Angular velocity + + ## Gyro offset + self.theta_offset = 0 self.dx = dx # Horizontal displacement of "cart" from center self.a_cart = [0] # Acceleration of cart @@ -49,9 +55,16 @@ class Pendulum: ## PID variables self.pid = True - self.kp = 1.3 + self.pidDelay = 0 + # self.kp = 1.3 + # self.ki = 0.0 + # self.kd = 0.1 + self.kp = 7.5 self.ki = 0.0 - self.kd = 0.1 + self.kd = 1.15 + + self.kp_m = 0.07 + self.kd_m = 0.1 def update(self, dt): self.doMath(dt) @@ -60,7 +73,12 @@ class Pendulum: ) if self.pid: - self.pidControl() + self.pidDelay += dt + + if self.pidDelay > 21: + self.pidControl(self.pidDelay) + self.theta_offset += self.pidDelay * (C_ANG_RES_LIMIT / 1000) ## One degree a second + self.pidDelay = 0 if abs(self.theta[self.index]) == C_FALL_ANG: self.fallen = True @@ -89,14 +107,14 @@ class Pendulum: def doMath(self, dt): ### ANGLE ### - ang_term1 = self.a_cart[self.index] * math.cos(self.theta[self.index]) - ang_term2 = self.v_cart[self.index] * math.sin(self.theta[self.index]) + ang_term1 = self.a_cart[self.index] * math.cos(self.theta[self.index] + self.theta_offset) + ang_term2 = self.v_cart[self.index] * math.sin(self.theta[self.index] + self.theta_offset) ang_term3 = ( self.v_cart[self.index] # Previous cart velocity * self.v_ang[self.index] # previous angle velocity - * math.sin(self.theta[self.index]) # Sin previous angle + * math.sin(self.theta[self.index] + self.theta_offset) # Sin previous angle ) - ang_term4 = C_GRAVITY * math.sin(self.theta[self.index]) + ang_term4 = C_GRAVITY * math.sin(self.theta[self.index] + self.theta_offset) # Angular acceleration self.a_ang.append( @@ -110,10 +128,13 @@ class Pendulum: ) # Angular displacement - self.theta.append( - self.theta[self.index] # Previous angle - + (self.v_ang[self.index + 1] * (dt / 1000)) - ) + theta_res = self.theta[self.index] + (self.v_ang[self.index + 1] * (dt / 1000)) + theta_round = theta_res % C_ANG_RES_LIMIT + if theta_round > 0.5 * C_ANG_RES_LIMIT: + theta_res = theta_res + theta_round + else: + theta_res = theta_res - theta_round + self.theta.append(theta_res) # Limit fall of pendulum self.theta[self.index + 1] = self.clamp( @@ -135,13 +156,20 @@ class Pendulum: ) # Cart acceleration - self.a_cart.append((-cart_term1 + cart_term2) / (2 * self.mass)) + a_res = (-cart_term1 + cart_term2) / (2 * self.mass) + self.a_cart.append(self.clamp(a_res, -C_A_CART_LIMIT, C_A_CART_LIMIT)) + # if abs(a_res) < C_A_CART_LIMIT: + # self.a_cart.append(a_res) + # else: + # self.a_cart.append(C_A_CART_LIMIT * (a_res / abs(a_res))) # Integrate acceleration to get velocity - self.v_cart.append( - self.v_cart[self.index] # Previous velocity - + (self.a_cart[self.index + 1] * (dt / 1000)) - ) + v_res = self.v_cart[self.index] + (self.a_cart[self.index + 1] * (dt / 1000)) # Previous velocity + + if abs(v_res) < C_V_CART_LIMIT: + self.v_cart.append(v_res) + else: + self.v_cart.append(C_V_CART_LIMIT * (v_res / abs(v_res))) # Cart displacement self.s_cart.append( @@ -167,6 +195,7 @@ class Pendulum: self.a_cart = [0] self.v_cart = [0] self.s_cart = [0] + self.theta_offset = 0 self.fallen = False self.update(0) @@ -193,7 +222,13 @@ class Pendulum: plt.show() - def pidControl(self): - error = self.theta[self.index] - result = (self.kp * error) + (self.ki * sum(self.theta)) + (self.kd * self.v_ang[self.index]) - self.a_cart[self.index] = result * 10 \ No newline at end of file + def pidControl(self, dt): + error = self.theta[self.index] + dt = (dt/1000) + # result = (self.kp * error) + (self.ki * sum(self.theta)) + (self.kd * self.v_ang[self.index]) #PID + result = (self.kp * error) \ + + (((error - self.theta[self.index - 1]) / (dt + 1 / 1000)) \ + * self.kd) \ + + (self.kp_m * self.s_cart[self.index]) \ + + (self.kd_m * self.v_cart[self.index]) + self.a_cart[self.index] = self.clamp(result * 10, -C_A_CART_LIMIT, C_A_CART_LIMIT) \ No newline at end of file diff --git a/src/sim/sim.py b/src/sim/sim.py index eabaaf4..0b5d1be 100644 --- a/src/sim/sim.py +++ b/src/sim/sim.py @@ -25,7 +25,7 @@ ui = SimUI(screen, pole) # Pendulum setup # Start angle in radians, length, mass, color -pendulum = Pendulum(0, 2, 0, 0.25, "red") +pendulum = Pendulum(0, 0.2, 0, 0.7, "red") pendulum.reset() dx = 0 # x offset dt = 1 # delta time @@ -100,9 +100,9 @@ while running: ui.grid(50, 0, 15) # Update PID values - pendulum.kp = slider_kp.getValue() - pendulum.ki = slider_ki.getValue() - pendulum.kd = slider_kd.getValue() + # pendulum.kp = slider_kp.getValue() + # pendulum.ki = slider_ki.getValue() + # pendulum.kd = slider_kd.getValue() # Update pendulum if not pendulum.fallen: