Sisällysluettelo

Avainsanoja: python, sympy, symbolinen laskenta, differentiaaliyhtälöt, numeerinen simulointi, fysiikka, mekaniikka, Lagrange mekaniikka, ongelmalähtöinen oppiminen

Vaikkei siltä ehkä näytä, tämä dokumentti on myös python-kielinen ohjelma, jonka alussa on hyvä kertoa allaolevalla tavalla, millaisen koodiston mukaan pythonin pitää tulkita dokumentin merkit.

In [1]:
# -*- coding: utf-8 -*-

Tutkiva oppiminen

Tavoitteeni oli pienellä demolla osoittaa,

  • miten tietotekniikalla voi tukea tutkivaa oppimista
  • että oppimistavoitteita kannattaisi päivittää vastaamaan tietotekniikan tarjoamia mahdollisuuksia

Heti projektin aluksi kävi ilmi, että joudun itsekin opiskelemaan lisää, joten tästä jutusta tuli kaksitasoinen: Välillä kirjoitan oppimisen ohjaajan näkökulmasta, välillä kuvaan omakohtaista tutkivaa oppimista. Koska opiskelin omin päin, opin kantapään kautta millaista apua opiskelija tarvitsisi.

Työelämässä uusia asioita joudutaan usein opettelemaan ja kehittämään omin päin — joskus yksin, onneksi useimmiten porukalla, mutta oppilaitoksissa tutkiva oppiminen edellyttää asiantuntevaa ohjausta. Oppimisen ohjaaja ei saa antaa valmiita ratkaisuja eikä mennä neuvomaan kuin tarvittaessa, mutta opiskelijoita ei saa jättää hakkaamaan päätään seinään. Se on liian hidas, tehoton ja turhauttava tapa opiskella.

Ohjaajan pitää osata opastaa, miten selvitä hankaluuksista ja mistä ja miten hakea lisätietoa. Valmiiden ratkaisujen tarjoaminen ei opeta sitä kaikkein tärkeintä, omin päin selviämistä. "Oppimaan ohjaaminen" vaatii erilaista kokemusta, taitoa ja asennetta kuin "perinteinen opettaminen".

Nosturin kuva
Kuva 1: Kiskoilla kulkevaan vaunuun, jonka massa on $m_{car}$, on ripustettu vaijerilla massa $m_{load}$. Voima $F_{car}$ liikuttaa vaunua ja voima $F_{load}$ kelaa kuorman kannatinvaijeria sisään ja ulos. $L$ on vaijerin pituus ja $\theta$ vaijerin poikkeama pystyasennosta.

Ongelmalähtöinen oppiminen

Olen ymmärtänyt, että ilmiöpohjaisessa oppimisessa tutkitaan jotain ilmiötä: Ei ratkaista annettua tehtävää vaan asetetaan itse kysymykset, joihin tuntuu tärkeältä löytää vastaus.

Minun lähtökohtani oli toinen, tavallaan vaatimattomampi. Halusin näyttää, miten matematiikkaa ja sen soveltamista, mekaniikkaa ja ohjelmointia voi oppia integroidusti ratkaisemalla sopivia tehtäviä. Ehkä voisin väittää, että tutkin tietynlaisia mekaniikan ilmiöitä konkreettisen esimerkin avulla.

Halusin kokeilla, miltä tuntuu matematiikan, tietotekniikan ja mekaniikan ongelmalähtöinen oppiminen: Aletaan ratkaista jotain konkreettista ongelmaa ilman tarvittavia esitietoja ja opiskellaan uusia asioita sitä mukaa kuin niiden tarve tule vastaan ja siinä laajuudessa, kuin käsillä oleva tehtävä vaatii.

Konkreettisella ongelmalla tarkoitan tehtävää, jonka tavoite on ymmärrettävissä. Mekaniikka tarjoaa hyviä matematiikan ja ohjelmoinnin harjoitustehtäviä, koska ongelmat ovat niin konkreettisia, että opiskelija voi ymmärtää, mitä on ratkaisemassa. Silti yksinkertaisenkin mekaanisen järjestelmän toiminta voi olla yllättävän monimutkaista ja siksi mielenkiintoinen ja opettavainen tutkimuskohde. Hankkeeseen kuuluisi tietysti myös havainnekuvan piirtäminen järjestelmästä ja matemaattista tekstiä sisältävän dokumentin tuottaminen. Ellei piirtäminen onnistuisi, pysähdyttäisiin harjoittelemaan sitä.

Minulle videopelimäisen simulaattorin tekeminen ylläolevan kuvan nosturista oli riittävän konkreettinen ja mielenkiintoinen tehtävä. Arvelin, että opiskelijat jonkin aikaa nosturin kuljettajaa leikittyään kiinostuvat kehittämään algoritmia, joka ratkaisee, miten siirtää kuorma automaattisesti paikasta toiseen.

Yksi syy Angry Birdsin viehätykseen on kuulemma se, että olioiden ja kappaleiden liikkeet noudattavat fysiikan lakeja. Mielenkiinnon herättämiseksi opiskelijat voisi laittaa arvioimaan Angry Birdsiä tuosta näkökulmasta. Angry Birds -fysiikkaa

Itse haluaisin kokeilla, miten tämän nosturihankkeen myötä oppimillani taidoilla voisi tutkia vaikka laskettelurinteiden hyppyreiden muotoja ja optimaalista ponnistusta. Eli kaikki ehdotukseni harjoitustöistä olisivat kovin samanlaisia eivätkä kiinnostaisi pätkääkään isoa osaa opiskelijoista.

Sopivien aiheiden ideointi pitäisi tehdä yhdessä opiskelijoiden kanssa. Aluksi opettajan pitäisi selittää, millaisia taitoja harjoitustyön myötä on tarkoitus oppia. Oppimistavoitteet pitää ymmärtää väljästi — perinteinen yksityiskohtainen luettelo ei johda tutkivaan oppimiseen eikä käytännön kannalta oleellisten taitojen oppimiseen.

Eräs blogikirjoittaja selitti oppimista käyttäen jääkiekkoa esimerkkinä.

Jääkiekkovalmentajan pitää osata "lukea peliä" kyetäkseen opastamaan joukkuetta. Sääntöjen osaaminen ei riitä. Ei sekään, että osaa luistella. Pelaajakokemuksesta on hyötyä, mutta pelatessa peliä ei ehtine tarkastella valmentajan näkökulmasta. Pelien katsominen opettanee lukemaan peliä, mutta viime kädessä vasta valmentajana toimiminen opettaa hakemaan vastausta kysymykseen, miten auttaa oma joukkue voittoon.

Miltään kurssilta ei valmistuta huippuvalmentajaksi eikä pelaajaksi:

  • Jääkiekko
    • säännöt: 9op, arvosana 2
    • luistelu: 12op, arvosana 4
    • mailan käyttö: 9op, arvosana 3
    • syöttäminen: 5op., arvosana: 5
    • laukaisut maaliin 8op., arvosana 4
  • Yhteensä:

(5: NHL-tason osaaminen, 4: maajoukkeutason osaaminen, ...)

Esimerkki harjoitustehtävästä

Tehtävänannosta tulee selkeämpi, jos heti aluksi selittää, mikä on tavoite. Esimerkiksi tässä nosturiprojektissa tavoitteena on tutustua yksinkertaisen mekaanisen järjestelmän toimintaan ja menetelmiin, jolla sitä voi analysoida.

Yksityiskohtamisempi tehtävä on kehittää menetelmiä nosturin ohjaamisen avuksi. Tämän tavoitteen voi esitellä videon avulla. Oheisessa videossa siirtelen ensin nosturia käsiohjauksella, sen jälkeen kytken päälle säädön, joka pitää nostovaijerin vakiomittaisena. Lopuksi annan automaattiohjauksen siirtää kuorman paikasta toiseen.

Demovideo

Koska videolta ei näe kaikkea, pitää tietysti opetella tuottamaan esimerkiksi allaolevan kaltaisia graafisia esityksiä järjestelmän toiminnasta:

Minun näkemykseni tavoitteesta lienee tullut selväksi, vaan tulinko antaneeksi aiemmin tuomitsemani valmiin ratkaisun?

Nosturiharjoituksestani ei taida tulla ilmiöoppimista eikä oikein tutkivaa oppimistakaan vaan matematiikan, fysiikan ja ohjelmoinnin opiskelua integroiva tarkasti ohjeistettu harjoitustyö.

Ilmiöoppiminen voisi alkaa esimerkiksi videoimalla kaikkea, missä arvelee fysiikan lakien vaikuttavan: nostureita, autoja, kiven heittämistä, puun heilumista, meren aaltoja, pilviä, ...

Oman kokemukseni mukaan ohjelmointia ja matematiikan soveltamista voi oppia muokkaamalla muiden ohjelmia ja ratkaisuja itselleen sopivaksi. Niin tässä tein itsekin. Jos olisin aloittanut kaiken alkeista, ei tästä hommasta olisi tullut mitään.

Tämän harjoituksen myötä aloin ymmärtää, miksi tutkivasta oppimisesta tuleekin niin helposti vallan muuta. Ymmärrän, miksi opiskelijat niin helposti pitävät tutkivaa oppimista keinotekoisena. Pahimmassa tapauksessa omalle ajattelulle ei jää tilaa: "Tutkikaa tätä ongelmaa, tällä menetelmällä ja päätykää tähän tulokseen."

Ilmiöoppimminen, tutkiva oppiminen ja ongelmalähtöinen oppiminen ovat hyvä tapa oppia, mutta vaatii oppimisen ohjaajalta taitoa ja kokemusta. Helposti nurjahtaa suorittamaan perinteistä yksityiskohtaiset oppimistavoitteet luettelemalla koottua suunnitelmaa. Ne tuntuvat yhtä välttämättömiltä kuin tutkivan oppimisen keinoin opittavat aivan muunlaiset taidot.

Miltä tuntui, miten kävi

Tieteellisissä artikeleissa johtopäätökset esitellään lopuksi. En esitä selkeitä johtopäätöksiä, koska tämä ei ole sellainen jämäkkä tutkimus, jolla niitä voisi perustella. Sensijaan kirjaan tähän johdannon tapaiseksi kokemuksiani ja tuntemuksiani.

Varasin tähän projektiin viikon verran aikaa, mutta loppujen lopuksi se vei reilun kuukauden. Jään tällaisiin harjoituksiin koukkuun kuin jotkut videopeleihin. Osan ajasta työskentelin flow-tilassa ajan kulua huomaamatta, mutta stressiäkin koin: Aika tuntui kuluvan turhanaikaiseen puuhasteluun, mutta toisaalta kaiken aherruksen jälkeen halusin jotain viimeisteltyä esitettäväksi.

Halusin ratkaista ongelman aina vain tyylikkäämmin ja perusteellisemmin, vaikka toisaalta mielenkiintoni kaikkoaa, kun itse näen tuloksen, josta näkee, että ratkaisuni toimii periaatteessa.

Olen innoissani oivallettuani selityksen jollekin asialle ja haluan kertoa sen kaikille, mutta riittävän perusteellisen selityksen laatiminen on enimmäkseen tylsää rutiinia. Selitystä laatiessaan huomaa, kuinka moni "itsestäänselvyys" pitäisi selittää. Kunnollinen selitys vaatisi esimerkin, piirroksia, demonstraatiota, — —

Vaikka minulla oli hyvät perustiedot niin mekaniikasta, tarvittavasta matematiikasta kuin ohjelmoinnistakin, kaikkea jouduin opiskelemaan lisää.

Kokoilin googlaimella eri lähteistä "reseptin", jonka avulla sain tehtävän ratkaistua, vaikken aina kunnolla ymmärtänyt, mihin resepti perustui. Kun aloin kirjoittaa tätä juttua, palasin lähteisiin ja ymmärsin monta asiaa vähän paremmin.

Ohjelmoinnin osalta opettelin aiemmin osaamaani elegantimpia tapoja koodata algoritmit. Kirjoitin uudelleen toimiviakin osia koodista, koska ne alkoivat näyttää kömpelöiltä. Tällaisessa projektissa koodia joutuu muuttamaan jatkuvasti, joten koodista on pakkokin tehdä korkealaatuista, koska huonosta koodista joutuu kirjoittamaan uudelleen paljon suuremman osan kuin hyvin koodatusta ja huonon koodin muokkaaminen on vaikeaa ja johtaa vaikeasti löydettäviin virheisiin. Tällainen projekti on hyvä tilaisuus oppia ymmärtämään, millaista on hyvä koodi.

Joitain asioita olen selittänyt hyvinkin seikkaperäisesti, mutta jotkut vaikeatkin kohdat ovat vielä selitystä vailla. Ainakin toistaiseksi tämä juttu edellyttää esitietoina jonkin moisen ymmärryksen siitä, mitä ovat funktio ja sen derivaatta ja integraali. Ohjelmoinnistakin olisi hyvä olla pikkuisen kokemusta.

Kirjoitin tämän jutun tarpeisiin sovitetun selityksen funktioista ja differentiaaliyhtälöistä. Ei sekään yksistään riittäne pohjatiedoiksi. Kaiken kaikkiaan juttuihini pätee, että vaikeaselkoinen kohta on huonosti tai vajavaisesti selitetty. Virheitäkin olen saattanut tehdä. Virheitä ja huonosti selitettyjä kohtia on oppikirjoissakin, vaikka ne on tarkistettu paljon perusteellisemmin kuin tämä juttu.

Systemaattista kirjallisuushakua en tehnyt, mutta silmiini sattui mm. systemaattinen esittely ohjelmoinnista ja simuloinnista

Ongelmalähtöisessä oppimisessa pitää löytää toimiva perusteellisuuden taso. Kaikessa "insinöörilaskennassa" tulee vastaan muun muassa pythagoraan lause. Tämän harjoituksen kannalta riittää piirtää kuva suorakulmiosta ja kertoa, että $a² = b² + c²$. Matemaattisesta näkökulmasta pitänee opettaa myös termit hypotenuusa ja kateetti ja esittää lauseen todistus.

Jaarittelusta

Joutava jaarittelu rasittaa lukijaa, mutta liika tiiviys rajaa lukijat vain niihin, jotka jo tuntevat asian. Hyvää tekstiä jaksaa lukea pitkäänkin, huonon heittää sikseen nopeasti. Minä haluaisin selittää asiat ymmärrettävästi, vaikka pitemmästikin, mutta aika ja tarmo tuppasivat loppumaan kesken.

Joskus minua houkuttaa selittää asiat "epämatemaattisin" esimerkein ja selityksin, vaikka sellainen minua epäilyttääkin: Matematiikka on matematiikkaa, joten olisi kai hyvä osata selittää sitä matemaattisesti — johdattaa lukija matemaattiseen ajatteluun.

Matematiikka kokonaisuudessaan on järkeen käypää, kaikelle on ymmärrettävä täsmällinen perustelu. Matemaattisen tekstin lukeminen vaatii kuitenkin paljon harjaannusta. Minua on auttanut tulosten soveltaminen "käytäntöön". Ymmärrän kerta kerralta enenmmän lukemastani, kun tunnen hyvin jonkun ongelman, johon sitä voi soveltaa.

Matematiikka on eksakti tiede, jossa on tärkeää tehdä kaikki täsmälleen oikein. Tietyt matemaattiset totuudet pätevät vain tarkalleen tiettyjen ehtojen vallitessa. Kun tuntee hyvin matematiikkaa ja fysiikan lakeja, ymmärtää paremmin, miten erilaisia ongelmia voi ratkaista. Muutaman kerran oivalsin vasta parin päivän puurtamisen jälkeen, että yritin ratkaista ongelmaa aivan väärällä tavalla. Matemaattisemmin ajatteleva olisi varmasti heti huomannut virheeni samalla lailla kuin fysiikkaa ymmärtävä tunnistaa yrityksen tehdä ikiliikkuja. Ikiliikkujaa kannattaa yrittää rakentaa vasta, jos on ensin löytänyt virheen nykyfysiikasta.

Tutkiva oppiminen vs. ainelähtöinen oppiminen

Koulun fysiikan tunneilla moni olennainen ilmiö jäi selittämättä tai selitettiin tosi oudosti, koska matematiikan tunneilla ei (vielä) oltu ehditty käsitellä tarvittavia asioita. Differentiaalilaskentaa ei voitu hyödyntää ennen kuin sitä oli käsitelty matemaattisen perusteellisesti matematiikan tunneilla. (Kävin koulua 60 ja 70-luvuilla)

Osaaminen soveltaa eri aineiden tietoja ja taitoja ei kehity kuin harjoittelemalla. Jos "sovellusesimerkit" räätälöidään ainekohtaisesti yksittäisen opitun taidon harjoitteluun ne ovat väistämättä rajoittuneita ja keinotekoisen tuntuisia.

Vaatimus tutkia pelkästään aitoja elävän elämän ilmiöitä voi puolestaan johtaa liian isoihin projekteihin, mikäli niistä halutaan jotain oppiakin. Tärkeintä on, että tutkittavat ilmiöt ovat oikealla tavalla opettavaisia. (Mitä tuo tarkoittaa, vaatii huolellista miettimistä.)

Aloin kirjoittaa tätä juttua sillä ajatuksella, että lukijalle melkein kaikki saattaa olla uutta, joten selitän kaiken tarvittavan ainakin pintapuolisesti. Pian minulle valkeni, että kirjoitettavaa tulisi satoja sivuja monista eri aiheista vaikka vain pintaa raapaisisi. Wikipediasta löytyy artikkeli kaikesta tarvittavasta, mutta niihin viittaaminen ei riitä — ne eivät sovi itseopiskeluun. Tämänkin takia ilmiöoppiminen ja ongelmalähtöinen oppiminen edellyttävät osaavaa ohjausta.

Korkeakoulupiskelua aloittaessani meille sanottiin, että viiden vuoden kuluttua ymmärrätte, mitä hyötyä opiskeltavista aineista on. Viisi vuotta on pitkä aika. Oppien soveltamista pääsi pienimuotoisesti harjoittelemaan laboratoritöissä, mutta laajemmin vasta diplomityötä tehdessään. Aikaansaamista ja teoreettisten tietojen ja taitojen soveltamista olisi hyvä päästä harjoittelemaan jo opintojen varhaisessa vaiheessa sekä motivaation että ammattitaidon kehittymisen kannalta.

Soveltavien monialaisten harjoitustöiden tekeminen kestää kauan, kun samalla opiskellaan tarvittavia "esitietoja". Silti se tuntuu motivoivammalta kuin "esitietojen" tahkoaminen vuositolkulla ennen kuin pääsee näkemään, mihin niitä voi soveltaa. Ehkä "esitietoihin" kannattaisi syventyä perusteellisemmin vasta muutaman soveltavan harjoituksen jälkeen: Tehdään ensin ohjaajan avulla, sitten pohditaan, mitä oikein tuli tehtyä ja lopuksi selvitetään teoreettiset taustat ja perustelut.

Ryhmissä tekemällä oppiminen

Yhdessä tekeminen on tärkeä taito, jota ei voi oppia kuin tekemällä hommia ryhmissä. Lisäksi, jos ryhmä toimii hyvin, se jaksaa yrittää pidempään kuin yksinäinen puurtaja ja pärjää pidemmälle omin voimin ja vähemmällä opastuksella.

Kun opiskelu perustuu tekemiseen, yksilöt ja ryhmät etenevät omaa tahtiaan ja vain harvoin kaikki tarvitsevat samanlaista ohjausta. Kullakin on oman laisensa ongelmat tai ainakaan sama selitys ei auta kaikkia. Koko porukalla luennoinnista ei siksi ole paljon hyötyä, vaan oppimisen ohjaaja joutuu antamaan henkilö- ja ryhmäkohtaista ohjausta. Koska oppimisen ohjaajan aika on rajallista, kaikille tarpeellinen yleinen ohjeistus on oltava webissä kaikkien saatavilla juuri silloin, kun he sitä tarvitsevat. Yleinen ohjeistus ei saisi olla valmiin ratkaisun antamista, mikä tekee ohjeiden kirjoittamisesta haastavaa. Se kun vie muutenkin valtavasti aikaa.

Ylläoleva sopii melko tiukasti rajattuihin harjoituksiin. Varsinaisessa ilmiöoppimisessa opiskelijat voivat suuntautua enemmän oman mielenkiintonsa mukaan, jolloin tehtävän tekemiseen ei voine edes olla valmiita ohjeita.

Ilmiöoppiminen ja tutkiva oppiminen tuntuvat minusta juuri oikeilta tavoilta oppia, mutta mitä tarkemmin asiaa miettii, sitä paremmin huomaa, miten vähän oikeasti ymmärtää, kuinka tarkkana pitää olla siinä, miten tällaisen opetuksen järjestää ja miten sitä ohjaa.

Tietokoneet oppimisen apuna

Pikaisella googlailulla löysin monta esimerkkiä, jotka auttoivat minua ratkaisemaan nosturin ohjaustehtävän. Näytti kuitenkin siltä, että kukaan ei ole innostunut opetustarkoituksessa simuloimaan esimerkkeinä käyttämiään järjestelmiä. Vaikutti siltä, että opetuksessa käytettävien esimerkkien tehtävät oli keksitty sopivasti niin, että ne saattoi ratkaista analyyttisesti, ei niin, että tulos olisi ollut — ainakaan minun mielestäni — mielenkiintoinen. Lisäksi vaikutti siltä, että samoja esimerkkejä kierrätetään uudistamatta niitä vastaamaan nykyisten tietoteknisten työkalujen tarjoamia mahdollisuuksia.

Minun kouluaikoinani dynaamisten järjestelmien simulointi oli kaiketi mahdollista vain hyvin varustetuissa yliopistoissa ja tutkimuslaitoksissa, nyt saman voi tehdä läppärillä. Miksei siis jo lukiossa simuloitaisi järjestelmien toimintaa? Fysiikan tunneilla harjoiteltiin voimien laskemista, muttei voitu tutkia, millaista liikettä voimat saavat aikaan vähänkään monimutkaisemmassa järjestelmässä. Mielenkiintoisin osuus jäi näkemättä. Ainakin minun mielestäni mielenkiintoisin.

Matematiikan kursseilla tarkasteltiin differentiaaliyhtälöistä vain analyyttisesti käsityönä ratkaistavissa olevia, mikä rajoitti sovellusesimerkit kaikkein yksinkertaisimpiin. Matematiikassa tärkeintä tietysti onkin oppia ymmärtämään differentiaalilaskennan perusteita, mutta minusta mielenkiintoisinta on päästä soveltamaan matematiikka ja nähdä, miten esimerkiksi fysiikan lait määräävät sen, mitä oikeasti tapahtuu.

En tiedä, kuuluvatko differentiaaliyhtälöt nykyiseen matematiikan opetussuunnitelmaan. Ellei, siihen lienee hyvä syy. Jotkut lienevät miettineet, miten parhaiten oppia itse matematiikkaa. Yksittäisten laskumenetelmien hallitseminen ei välttämättä ole siinä tärkeintä. Kun hallitsee hyvin itse matematiikan perusteet, erilaisten matemaattisten menetelmien oppiminen ei ole iso vaiva.

Fysiikka on vaikea ymmärtää, ellei ymmärrä differentiaaliyhtälöitä, derivointia ja integrointia. Ehkä tarvittaviin matemaattisiin työkaluihin voitaisiin tutustua tarpeen mukaan.

Tässä jutussa esittelen "insinöörilaskentoa", mutta sitä ei olisi olemassa ilman matematiikkaa. Insinöörilaskennalla tarkoitan jonkun johtamien matemaattisten tulosten soveltamista kuin "resepteinä" perehtymättä siihen, miten tulokset on johdettu ja mihin ne perustuvat. Näin pääsee nopeasti kokeilemaan mielenkiintoisten ongelmien ratkaisemista, mutta voi myös erehtyä yrittämään mahdottomia tai epäonnistua siksi, että on valinnut väärän "reseptin". Tosielämässä pitää perehtyä "reseptien" taustoihin sen varmistamiseksi, että ratkaisu pätevä kaikissa mahdollisissa tilanteissa, missä reseptiä soveltaa.

Projektissa tarvittavat ohjelmat

Jos haluat tehdä itse jonkin vastaavan harjoituksen tai kokeilla esittelemiäni menetelmiä, tarvitset python-ohjelmointiympäristön. Pythonista on olemassa versiot python2 ja python3, jotka eivät ole aivan yhteensopivia. Käytän python3:a.

Opettelin ohjelmoimaan pythonilla pari vuotta sitten, joten kaikkea en vielä osaa tehdä järin tyylikkäästi.

Asensin pythonin tietokoneeseeni linuxin pakettienhallintajärjestelmällä. Käytän spyder-nimistä ohjelmointiympäristöä, mutta muitakin hyviä on. Pythonin voi tietysti asentaa myös windowsiin ja appleen.

Pythonin asennuksen myötä koneelleni ilmestyivät komennot pip, pip2 ja pip3 pip3.6. Ne ovat pythonin lisämodulien hallintaohjelmia. Niillä voit installoida sellaisia python-moduleita, joita ei löydy koneesi pakettienhallinnasta. Tässä harjoituksessa tarvittavia lisämoduleita ovat esimerkiksi sympy ja numpy. (Jos installoit numpy-paketin pakettienhallintajärjestelmällä, installoi myös numpy-devel.)

Numeerisessa optimoinnissa käyttämäni scikits.bvp_solver ohjelmistopaketin asensin komennolla pip3.6 install scikits.bvp_solver, koska sitä ei löytynyt linux-järjestelmäni pakettienhallinnasta.

(Harjoituksen vuoksi latasin osoitteesta https://pypi.python.org/pypi/scikits.bvp_solver lähdekielisen version. Komennoilla python3.6 setup.py build ja python3.6 setup.py install sain asennettua python3-version scikits.bvp_solversta.)

Mikäli tämän dokumentin nimen loppuliite on html luet html-muotoon tallennettua jupyter notebookia. Jos haluat tutustua alkuperäiseen jupyter notebookiin installoi jupyter paketti joko linuxin pakettienhallinnalla tai pip3 install jupyter -komennolla. Kun tämän jälkeen annat komennon jupyter3-notebook, selaimeesi aukeaa notebook sovellus. (notebook tiedosto ei aukea selaimeen edelläolevaa linkkiä klikkaamalla, mutta voit downloadata sen klikkaamalla hiiren oikeaa.)

Notebookissa on kahdenlaisia soluja: code ja Markdown. code-soluissa on suoritettavaa python-koodia, Markdown-soluissa tämänkaltaista dokumentoivaa tekstiä. Tuplaklikkaamalla Markdown-solua näet, missä muodossa teksti on kirjoitettu ja pääset editoimaan sitä.

Jos itse alat kirjoittaa jupyter-notebook tekstiä, kannattaa tutustua Markdown notaatioon.

Jos haluat kirjoittaa matemaattisia lausekkeita, kannattaa tutustua LaTeX komentoihin. Myöhemmin tässä dokumentissa näet niistä esimerkkejä.

Jupyter-notebookiin tutustuin tätä projektia aloitellessani ja olen opetellut käyttämään sitä tätä hommaa tehdessäni. Opin koko ajan lisää hienouksia, joita tulee kiusaus kokeilla ja esitellä. Loputon suo sekin :-)

zip-arkisto kaikesta tämän projektin koodista

Siltanosturi — simulaattori ja liikeradan optimointi Pythonilla

Alustuksia

Python-ohjelmien alussa luetellaan python-moduleja, joiden funktiota käytetään hyväksi.

Pyörittelen tässä ohjelmassa symbolisia lausekkeita, joista generoidaan python-koodia numeerista laskentaa varten. Koodista tulee niin monimutkaista, etteivät sen lausekkeet mahdu yhdelle riville. textwrap-paketista löytyvä metodi fill osaa katkoa rivit automaattisesti pythonin hyväksymällä tavalla.

... import fill tarkoittaa, että voin kutsua funktiota fill suoraan.

In [2]:
from textwrap import fill, wrap

from IPython.display import display

Numpy-moduli tarjoaa numeerisen laskennan apufunktioita. import numpy muotoinen import komento edellyttää, että esimerkiksi numpy-modulin funktioita array kutsuessani, minun pitää numpy.array

Tyylikkäintä olisi ollut käyttää vain yhtä import-muotoa, mutta tulin kokeilleeksi kaikkia vaihtoehtoja.

In [3]:
import numpy

import math # matemaattisia funktioita kuten sin(), cos()

Re-modulissa on regular expressions funktioita merkkijonojen käsittelyyn. Regular expressions on käytettävissä kaikissa kunnollisissa editoreissa ja kaikki ohjelmointikielet tukevat niitä. Regular expressions-funktioita on vaikea oppia käyttämään, mutta ne nopeuttavat editointia ja ohjelmoinnissa ei aina edes selviä ilman niitä.

In [4]:
import re

Symboliseen laskentaan käytetään sympy-modulia. Se ei kuulu pythonin perusasennukseen, joten installoin sen erikseen komennolla pip3 install sympy.

In [5]:
import sympy as sp

... as sp tarkoittaa, että voin tästä ohjelmasta kutsua sympy:n funktioita lisäämällä funktion nimen eteen sp. Tämä tekee koodista selkeämpää.

In [6]:
sp.init_printing()

Ylläoleva varmistaa, että matemaattiset lausekkeet tulostetaan niin hyvin, kuin tietokoneesi sallii. Olen installoinut LaTeX ohjelmiston yms. joten minulle yhtälöt tulostuvat siististi.

(init_printing()-funktiolle voi antaa monenlaisia parametreja, mutta niihin perehtyminen ei tuntunut hyödylliseltä.)

Symbolisten muuttujien esittely

Ohjelmoinnissa muuttujilla pitää antaa jokin arvo ennen kuin käyttää niitä.

In [7]:
massa = 3.0
nopeus = 10.0
KE = (massa * nopeus**2) / 2
print("kineettinen energia: ", KE)
kineettinen energia:  150.0

Symbolisessa laskennassa käytettävät muuttujat pitää määritellä symboleiksi. W_{kin} ja display(sp.Eq ovat temppuilua, jota en malttanut olla kokeilematta. Niitä ei ole tärkeä ymmärtää tässä vaiheessa.

In [8]:
massa, nopeus, ke = sp.symbols("m, v, W_{kin}")
KE = (massa * nopeus**2) / 2
print("kineettinen energia: ", KE)
# kumpikin seuraavista näyttää lausekkeen siistissä muodossa
display(sp.Eq(ke, KE))  # Temppuilua, unohda jos näyttää vaikealta
KE
kineettinen energia:  m*v**2/2
$$W_{kin} = \frac{m v^{2}}{2}$$
Out[8]:
$$\frac{m v^{2}}{2}$$

Ylimmällä rivillä yhtäsuuruusmerkin vasemmalla on siis python-koodissa käytettävät muuttujien nimet ja oikealla symbols -määrittelyssä se, miten ne näkyvät tulostettuina.

Useimmiten on tietysti selkeintä, että molemmat ovat samoja:

In [9]:
m, v = sp.symbols('m, v')
KE = (m*v**2)/2
KE
Out[9]:
$$\frac{m v^{2}}{2}$$

Joskus on hyödyllistä esitellä määrittelemätön funktio.

In [10]:
t = sp.symbols('t')
y = sp.Function('y')(t)
sol = sp.dsolve(sp.diff(y, t) - y)
display(sp.Eq(sp.diff(y, t), y))
display(sol)
$$\frac{d}{d t} y{\left (t \right )} = y{\left (t \right )}$$
$$y{\left (t \right )} = C_{1} e^{t}$$

Käskyllä y = sp.Function('y')(t) muuttujalle y annetaan arvoksi määrittelemätön funktio y(t).

Käskyllä sp.dsolve(sp.diff(y, t) - y) ratkaistaan differentiaaliyhtälö $\frac{d}{dt}y(t) = y(t)$

Muotoa sp.dsolve(sp.diff(y, t) = y) ei voi käyttää vaan on käytettävä joko muotoa sp.dsolve(sp.diff(y, t) - y) tai sp.dsolve(sp.Eq(sp.diff(y, t), y))

Yhtälöä a = b ei voi tulostaa sellaiseen display-funktiolla, vaan a = b pitää korvata lausekkeella sp.Eq(a,b).

Selityksen siihen, miksi yhtälöitä ei voi käyttää kuin esim. Macsymassa, olen lukenut, mutten kykene sitä tähän toistamaan.

Yhtälöryhmät, matriisit ja vektorit

(Tämä luku on vielä melkoinen huitaisu, mutta koska tässä projektissa pyöritellään matriiseja, on ne pakko edes mainita.)

Yksi yhtälö riittää, jos tarkastelemme esimerkiksi auton kiihdyttelyä ja jarruttelua suoralla tiellä, mutta nosturin kuorma liikkuu tasossa — oikealle ja vasemmalle sekä ylös ja alas.

Jos haluamme tarkastella kappaleen liikettä tasossa, meidän pitää laskea voimat ja kiihtyvyydet sekä x- että y-suunnassa. Yhtälöitä tulee kaksin kappalein, mutta ne ovat samanlaisia. Tällaisissa tapauksessa yhtälöiden kirjoittamista selkeyttävät matriisit ja vektorit.

Matriisi, jossa on kaksi riviä ja kaksi saraketta:

$ \underline A = \begin{bmatrix} a_{11} \ a_{12}\\ a_{21} \ a_{22}\\ \end{bmatrix}$

Matriisimuuttujan, kuten yllä A nimi on tapana alleviivata, $ \underline A$, tai kirjoittaa se lihavoituna, $ \mathbf{A}$. (Noudatan tätä tapaa aina kuin muistan ja jaksan.)

Vektori on matriisin erikoistapaus, jossa on vain yksi rivi tai sarake. Kaksi vektoria:

$\underline x = \begin{bmatrix} x_{1} \\ x_{2} \\ \end{bmatrix}$ ja $\underline b = \begin{bmatrix} b_{1} & b_{2} \end{bmatrix}$

Kahden matriisin yhtäsuuruus tarkoittaa sitä, että niissä vastaavissa kohdissa olevat alkiot ovat samansuuruisia.

Matriisien avulla voimme esittää lineaarisen yhtälöryhmän seuraavasti:

$$\begin{bmatrix} y_{1} \\ y_{2} \\ \end{bmatrix} = \underline y = \underline A \underline x = \begin{bmatrix} a_{11}x_{1} \ a_{12}x_{2}\\ a_{21}x_{1} \ a_{22}x_{2}\\ \end{bmatrix}$$

Kun minulle ensimmäisenä opiskeluvuonna alettiin tyrkyttää oppia matriiseista, pidin koko touhua aivan joutavana. Seuraava vuonna matriiseja kummitteli melkein jokaisen oppikirjan sivuilla.

Tässäkin projektissa esitän melkein kaikki yhtälöt matriiseina. Ohjelmointiakin helpottaa usein matriisien käyttäminen.

Matriisien käsittely pythonilla

Joissain tapauksissa matriisit voi esittää python-kielen listoina.

Pythonin listat pitää jollain tavalla esitellä/alustaa ennen niiden käyttöä. Turhan näköisten komentojen ainoa tarkoitus on kertoa listan pituus. Sama pätee myöhemmin esiteltäville matriiseille.

Python-kielessä indeksointi alkaa 0:sta, eli listan ensimmäisen alkion indeksi on nolla.

In [11]:
lista = [1,2,3,4]  # Esitellään neljän alkion mittainen lista 
print(lista)
print(lista[0])
print(lista[3])
[1, 2, 3, 4]
1
4
In [12]:
# listan voi luoda näinkin
lista = [i for i in range(1,5)]
lista
Out[12]:
$$\left [ 1, \quad 2, \quad 3, \quad 4\right ]$$

Kaksiulotteinen matriisi listojen listana. Laitan alkioiksi merkkijonoja, joista näkyvät alkioiden indeksit

In [13]:
A = [["a"+str(i)+str(j) for j in range(0,3)] for i in range(0,3)]
print("matriisi: ",A)
print("alkio (1,1): ",A[1][1])
print("rivi 2: ", A[2])
matriisi:  [['a00', 'a01', 'a02'], ['a10', 'a11', 'a12'], ['a20', 'a21', 'a22']]
alkio (1,1):  a11
rivi 2:  ['a20', 'a21', 'a22']

Pythonin Sympy-paketin symbolisen laskennan funktiot vaativat matriisit esitettäväksi sympyn Matrix-muodossa.

In [14]:
B = sp.Matrix(A)
B
Out[14]:
$$\left[\begin{matrix}a_{00} & a_{01} & a_{02}\\a_{10} & a_{11} & a_{12}\\a_{20} & a_{21} & a_{22}\end{matrix}\right]$$
In [15]:
# sympy-matriisista voi kätevästi poimia alkion, sarakkeen tai rivin
display(B, B[1, 1], B[:, 1], B[1, :])
$$\left[\begin{matrix}a_{00} & a_{01} & a_{02}\\a_{10} & a_{11} & a_{12}\\a_{20} & a_{21} & a_{22}\end{matrix}\right]$$
$$a_{11}$$
$$\left[\begin{matrix}a_{01}\\a_{11}\\a_{21}\end{matrix}\right]$$
$$\left[\begin{matrix}a_{10} & a_{11} & a_{12}\end{matrix}\right]$$

Operaatiota matriiseilla:

Transpoosi, kertolasku ja summa

In [16]:
u = sp.Matrix([sp.symbols("u_" + str(i)) for i in range(1, 4)])
v = sp.Matrix([sp.symbols("v_" + str(i)) for i in range(1, 4)])
print("u ja v")
display(u, v)

ut = u.T
print("u:n transpoosi")
display(ut)

uvt = u * v.T
print("uvt: vektorien tulo")
display(uvt)

utv = u.T * v
print("utv: vektorien tulo")
display(utv)

udotv = u.dot(v)
print("udottv: vektorien pistetulo")
display(udotv)

btb = B * B
print("btb: matriisi kerrottuna itsellään")
display(btb)

bpb = B + B
print("bpb: matriisi lisättynä itseensä")
display(bpb)
u ja v
$$\left[\begin{matrix}u_{1}\\u_{2}\\u_{3}\end{matrix}\right]$$
$$\left[\begin{matrix}v_{1}\\v_{2}\\v_{3}\end{matrix}\right]$$
u:n transpoosi
$$\left[\begin{matrix}u_{1} & u_{2} & u_{3}\end{matrix}\right]$$
uvt: vektorien tulo
$$\left[\begin{matrix}u_{1} v_{1} & u_{1} v_{2} & u_{1} v_{3}\\u_{2} v_{1} & u_{2} v_{2} & u_{2} v_{3}\\u_{3} v_{1} & u_{3} v_{2} & u_{3} v_{3}\end{matrix}\right]$$
utv: vektorien tulo
$$\left[\begin{matrix}u_{1} v_{1} + u_{2} v_{2} + u_{3} v_{3}\end{matrix}\right]$$
udottv: vektorien pistetulo
$$u_{1} v_{1} + u_{2} v_{2} + u_{3} v_{3}$$
btb: matriisi kerrottuna itsellään
$$\left[\begin{matrix}a_{00}^{2} + a_{01} a_{10} + a_{02} a_{20} & a_{00} a_{01} + a_{01} a_{11} + a_{02} a_{21} & a_{00} a_{02} + a_{01} a_{12} + a_{02} a_{22}\\a_{00} a_{10} + a_{10} a_{11} + a_{12} a_{20} & a_{01} a_{10} + a_{11}^{2} + a_{12} a_{21} & a_{02} a_{10} + a_{11} a_{12} + a_{12} a_{22}\\a_{00} a_{20} + a_{10} a_{21} + a_{20} a_{22} & a_{01} a_{20} + a_{11} a_{21} + a_{21} a_{22} & a_{02} a_{20} + a_{12} a_{21} + a_{22}^{2}\end{matrix}\right]$$
bpb: matriisi lisättynä itseensä
$$\left[\begin{matrix}2 a_{00} & 2 a_{01} & 2 a_{02}\\2 a_{10} & 2 a_{11} & 2 a_{12}\\2 a_{20} & 2 a_{21} & 2 a_{22}\end{matrix}\right]$$

Muodostin yllä matriisit u = {{u}} ja v = {{v}} kimurantilla, mutta pythonmaisella tavalla. sp.symbols palauttaa palauttaa symbolisen muuttujan, jonka nimeksi tulee sen argumenttina oleva merkkijono. Pythonissa listoja käsitellään usein noin, että for-silmukka on "listan sisällä". Saman voisi tehdä tavallisessa for-silmukassa, mutta ylläoleva lienee laskennallisesti tehokkaampi ja kätevä, kun siihen tottuu.

Matriisin tai vektorin transpoosissa u.T {{ut}} rivit on käännetty sarakkeiksi tai toisinpäin.

Matriiseja voi kertoa tai laskea yhteen, kunhan dimensiot täsmäävät.

u.T*v = {{utv}}

u.dot(v) = {{udotv}}

u*v.T = {{uvt}}

B*B = {{btb}}

B+B = {{bpb}}

Matriisien derivaattoja

Jatkossa tarvitsemme monen muuttujan funktion gradienttia eli matriisia, jossa kukin alkio on funktion osittaisderivaatta yhden muuttujan suhteen. Lisäksi saatetaan tarvita jacobian matriisia ja hessian matriisia, joista näytän myöhemmin esimerkin.

On pitkä juttu selittää, mihin näitä tarvitaan. Vielä pitempi jutusta tulee, jos selittää kunnolla. Tyydyn pariin yleiseen toteamukseen, joiden jälkeen näytän niistä esimerkin

Numeeristen menetelmien "käyttöohjeissa" pyydetään kirjoittamaan funktiot, jotka palauttavat gradientin, jacobin matriisin tai hessian matriisin. Helpottaa, jos tietää, mitä ne tarkoittavat.

Funktion $y = f(x_i), i = 0..n$ gradientti tietyssä pisteessä kertoo, kuinka jyrkästi $y$ muuttuu kutakin muuttujaa $x_i$ muutettaessa. Kun gradientti on jossain pisteessä nolla, piste on yksi funktion ääripisteistä.

Numeerisin menetelmin haetaan usein funktion minimiä tai maksimia. Gradienttia käytetään apuna päätettäessä, mistä suunnasta ja kuinka kaukaa ääriarvoa haetaan.

Jacobian matriisi on vektoriarvoiselle funktiolle $\underline y = \underline f(x_i), i = 0..n$ sama kuin gradientti yksiarvoiselle funktiolle.

Hessian matriisi on gradientin jacobin matriisi. Selviää parhaiten allaolevasta esimerkistä. Numeeriset optimointimenetelmät tarvitsevat tätä, mutta laskevat sen itse, ellei käyttäjä ohjelmoi sitä valmiiksi. $ y = H(\underline x_0)\underline x$ on funtkion lineaariapproksimaatio pisteessä $\underline x_0$, mihin perustuu sen hyödyllisyys numeerisissa menetelmissä.

Funktion $f(\mathbf x)$ gradienttia voi kutsua $f$:n derivaataksi vektorin $\mathbf x$ suhteen ja merkitä $grad\,f = \frac{d}{d\mathbf x}f(\mathbf x) = f_{\mathbf x}(\mathbf x)$

In [17]:
f = sp.symbols("f")
xx = [sp.symbols("x_" + str(i)) for i in range(0, 3)]
gradfx = [0 for i in range(0, 3)]

for i in range(0, 3):
    gradfx[i] = sp.diff(f(xx[0], xx[1], xx[2]), xx[i])

display(gradfx)
$$\left [ \frac{\partial}{\partial x_{0}} f{\left (x_{0},x_{1},x_{2} \right )}, \quad \frac{\partial}{\partial x_{1}} f{\left (x_{0},x_{1},x_{2} \right )}, \quad \frac{\partial}{\partial x_{2}} f{\left (x_{0},x_{1},x_{2} \right )}\right ]$$
In [18]:
# Esimerkki gradientista

f = xx[0]**2 + xx[0] * xx[1] + sp.sin(xx[0]) * sp.cos(xx[2]) + 3 * xx[2]
dfdx = [0 for i in range(0, 3)]

for i in range(0, 3):
    dfdx[i] = sp.diff(f, xx[i])
display('xx ', xx, 'f ', f, 'df/dx ', dfdx)
'xx '
$$\left [ x_{0}, \quad x_{1}, \quad x_{2}\right ]$$
'f '
$$x_{0}^{2} + x_{0} x_{1} + 3 x_{2} + \sin{\left (x_{0} \right )} \cos{\left (x_{2} \right )}$$
'df/dx '
$$\left [ 2 x_{0} + x_{1} + \cos{\left (x_{0} \right )} \cos{\left (x_{2} \right )}, \quad x_{0}, \quad - \sin{\left (x_{0} \right )} \sin{\left (x_{2} \right )} + 3\right ]$$
In [19]:
xx = sp.Matrix([sp.symbols("x_" + str(i)) for i in range(0, 3)])
ff = sp.Matrix([sp.symbols("f_" + str(i)) for i in range(0, 3)])
jacob_ff = sp.eye(3)
for i in range(0, 3):
    for j in range(0, 3):
        jacob_ff[i, j] = sp.diff(ff[i](xx[0], xx[1], xx[2]), xx[j])
jacob_ff
Out[19]:
$$\left[\begin{matrix}\frac{\partial}{\partial x_{0}} \operatorname{f_{0}}{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial}{\partial x_{1}} \operatorname{f_{0}}{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial}{\partial x_{2}} \operatorname{f_{0}}{\left (x_{0},x_{1},x_{2} \right )}\\\frac{\partial}{\partial x_{0}} \operatorname{f_{1}}{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial}{\partial x_{1}} \operatorname{f_{1}}{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial}{\partial x_{2}} \operatorname{f_{1}}{\left (x_{0},x_{1},x_{2} \right )}\\\frac{\partial}{\partial x_{0}} \operatorname{f_{2}}{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial}{\partial x_{1}} \operatorname{f_{2}}{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial}{\partial x_{2}} \operatorname{f_{2}}{\left (x_{0},x_{1},x_{2} \right )}\end{matrix}\right]$$
In [20]:
Hessian = sp.eye(3)  # Matriisi pitää jotenkin esitellä/alustaa ennen käyttöä
for i in range(0, 3):
    for j in range(0, 3):
        Hessian[i, j] = sp.diff(gradfx[i], xx[j])
display(xx, gradfx, Hessian)
$$\left[\begin{matrix}x_{0}\\x_{1}\\x_{2}\end{matrix}\right]$$
$$\left [ \frac{\partial}{\partial x_{0}} f{\left (x_{0},x_{1},x_{2} \right )}, \quad \frac{\partial}{\partial x_{1}} f{\left (x_{0},x_{1},x_{2} \right )}, \quad \frac{\partial}{\partial x_{2}} f{\left (x_{0},x_{1},x_{2} \right )}\right ]$$
$$\left[\begin{matrix}\frac{\partial^{2}}{\partial x_{0}^{2}} f{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial^{2}}{\partial x_{0}\partial x_{1}} f{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial^{2}}{\partial x_{0}\partial x_{2}} f{\left (x_{0},x_{1},x_{2} \right )}\\\frac{\partial^{2}}{\partial x_{0}\partial x_{1}} f{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial^{2}}{\partial x_{1}^{2}} f{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial^{2}}{\partial x_{1}\partial x_{2}} f{\left (x_{0},x_{1},x_{2} \right )}\\\frac{\partial^{2}}{\partial x_{0}\partial x_{2}} f{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial^{2}}{\partial x_{1}\partial x_{2}} f{\left (x_{0},x_{1},x_{2} \right )} & \frac{\partial^{2}}{\partial x_{2}^{2}} f{\left (x_{0},x_{1},x_{2} \right )}\end{matrix}\right]$$

Siltanosturin tilayhtälön johtaminen

Miten saada selkoa siitä, miten kuvan 1 siltanosturi reagoi sitä ohjaaviin voimiin?

Kuvittele itsesi nosturinkuljettaksi, jonka pitäisi siirtää nosturin koukussa riippuva kuorma rakenteilla olevan talon katolle. Mitä teet, riippuu ensinnäkin siitä, missä nosturin kuorma ja vaunu ovat. Ohjaamisen kannalta on oleellista myös, ovatko vaunu ja kuorma paikallaan vai liikkuvatko ne.

Nosturin ja vaunun paikka ovat nosturin tilasuureita. Tilayhtälöllä tarkoitetaan differentiaaliyhtälöä, joka kertoo, miten nosturin tila muuttuu ajan myötä ja miten ulkoiset ohjaukset vaikuttavat siihen.

Newtonin lain mukaan $F(t)=m\cdot a(t)$ missä $F(t)$ on kappaleeseen kohdistuva voima, $m$ kappaleen massa ja $a(t)$ kappaleen kiihtyvyys.

Kiihtyvyys on nopeuden muutosnopeus eli $a(t)=\frac{d}{dt}v(t)=\dot{v}(t)$ ja nopeus on paikan muutosnopeus eli $v(t)=\frac{d}{dt}x(t) = \dot{x}(t)$. Näistä saamme $a(t) = \frac{d^{2}}{dt^{2}}x(t) = \ddot{x}(t)$.

Voima $F(t)$ voi olla järjestelmän sisäinen voima, esimerkiksi heilurin painoa kannatteleva voima, tai ulkoinen voima, kuten esimerkiksi voimat $F_{car}$ ja $F_{load}$, joilla ohjataan kuvan nosturia.

Mekaanisten järjestelmien toimintaa siis kuvaa differentiaaliyhtälö

$$\frac{d^{2}}{dt^{2}}\mathbf{x}(t) = \mathbf{f}(\mathbf{x}(t),\mathbf{F}(t),t)$$

Yhtälö on usein usean muuttujan differentiaaliyhtälöryhmä. Siitä tarkemmin seuraavassa kappaleessa.

Jos haluaa nähdä, miten järjestelmä toimii, pitää ratkaista tuo differentiaaliyhtälö. Yksinkertaisessa tapauksessa differentiaalyhtälö voidaan ratkaista analyyttisesti eli sopivia matematiikan kaavojan soveltaen. Monimutkaisessa tapauksessa differentiaaliyhtälö täytyy ratkaista numeerisesti iteroimalla.

Teemme simulaattorin, jolla voimme leikkiä nosturinkuljettajaa. Tarvitsemme differentiaaliyhtäläiden numeerista ratkaisemista, koska

  • analyyttisessa ratkaisussa voimat pitää tietää etukäteen alkuhetkestä loppuhetkeen, mutta simulaattorissa ohjaaja muuttelee ohjauksia koko ajan
  • yhtälöt ovat liian monimutkaisia analyyttisesti ratkaistaviksi

Nosturin vaunun paikka voidaan ilmoittaa esimerkiksi sen etäisyytenä kiskojen vasemmasta päästä: $\text{ vaunun paikka } = x_{car}(t)$. Vaijerin kuorman paikan ilmoittamiseen tarvitaan kaksi muuttujaa — sen paikka vaakasuunnassa ja kuinka korkealla se on. Samoin kuorman nopeuden esittämiseen tarvitaan kaksi komponenttia. Paikan ja nopeuden voi esittää kätevästi vektroreilla:

$$\text{ kuorman paikka } = \begin{bmatrix} x_{load}(t)\\ y_{load}(t) \end{bmatrix}$$

$$\text{ kuorman nopeus } = \begin{bmatrix} v^{x}_{load}(t)\\ v^{y}_{load}(t) \end{bmatrix} = \begin{bmatrix} \frac{d}{dt}x_{load}(t)\\ \frac{d}{dt}y_{load}(t) \end{bmatrix}$$

$$\text{ ja kiihtyvyys} = \begin{bmatrix} a^{x}_{load}(t)\\ a^{y}_{load}(t) \end{bmatrix} = \begin{bmatrix} \frac{d}{dt}v^{x}_{load}(t)\\ \frac{d}{dt}v^{y}_{load}(t) \end{bmatrix}$$

Vektoriesitystä käytetään eräänlaisena lyhennysmerkintänä niin, että voimme käyttää samoja yhtälöitä riippumatta siitä, tarkastelemmeko 1,2 vai 3 -ulotteista avaruutta.

Vektorien avulla voimme merkitä $\underline{F}(t)=m \underline{a}(t) = m\underline{\dot{v}}(t) = m\underline{\ddot{x}}(t)$ niin tasossa kuin 2 tai 3 -ulotteisessakin avaruudessa.

Tästä saamme differentiaaliyhtälön $\underline{\ddot{x}}(t) = \underline{F}(t)/m$

Juonipaljastus: Homman edetessä saamme aikaiseksi differentiaaliyhtälöitä, joita ei voi ratkaista kuin numeerisilla menetelmillä. Differentiaaliyhtälöiden numeeriset ratkaisumenetelmät taas soveltuvat usein ainoastaan sellaisille differentiaaliyhtälöryhmille, joissa on vain ensimmäisiä derivaattoja. Niin minunkin valitsemani (= löytämäni).

Sijoituksilla $\underline{a}(t) \mapsto \underline{\dot{v}}(t) \text{ ja } \underline{v}(t) \mapsto \underline{\dot{x}}(t)$ saamme seuraavan matriisimuotoisen yhtälöryhmän:

$$\begin{equation} \begin{bmatrix} \frac{d}{dt}v_{car}(t)\\ \frac{d}{dt}v^{x}_{load}(t)\\ \frac{d}{dt}v^{y}_{load}(t)\\ \frac{d}{dt}x_{car}(t)\\ \frac{d}{dt}x_{load}(t)\\ \frac{d}{dt}y_{load}(t)\\ \end{bmatrix} = \begin{bmatrix} (F_{car}(t)-F^{x}_{load}(t))/M_{car}\\ F^{x}_{load}(t)/M_{load}\\ F^{y}_{load}(t)/M_{load}\\ v_{car}(t)\\ v^{x}_{load}(t)\\ v^{y}_{load}(t) \end{bmatrix} \end{equation} $$

Yhtälö kuvaa, miten nosturi reagoi sitä ohjaaviin voimiin $F_{car}$ ja $F_{load}$.

Lagrangen mekaniikka

Ylläolevassa yhtälössä esiintyvien voiman $F_{load}$ komponenttien $F^{x}_{load}(t) $ ja $F^{y}_{load}(t) $ ratkaiseminen näytti helpolta, mutta osoittautui minulle liian vaikeaksi. Merkit ehkä sain lopulta oikein, mutta kuinka suuri on vaijeria kuormittava voima missäkin tilanteessa? Sen sijaan, että olisin saman tien päässyt kirjoittamaan koodia jouduinkin opiskelemaan aivan uusia asioita mekaniikasta.

Googlaimella minulle selvisi, että Lagrangen mekaniikka saattaisi auttaa. Lagrangen mekaniikkaa sovellettaessa voimien lausekkeiden sijasta kirjoitetaan liike- ja potentiaalienergian lausekkeet, mikä sopivasti valitussa koordinaatistossa on paljon helpompaa kuin yrittää ratkaista voimien lausekkeita.

Myöhemmin opin, että järjestelmän sisäisiä voimia — vaijeria jännittävää voimaa — ei tarvitse miettiä lainkaan, mutta ulkoiset voimat — ohjausvoimat ja kitkavoimat — aiheuttivat kuitenkin pikkuisen päänvaivaa.

Lagrangen mekaniikan mukaan järjestelmän liike- ja potentiaalienergian erotukselle $L$ ja q-koordinaatistossa esitetyille ulkoisille voimille $Q$ pätee:

$$ \frac{d}{dt}(\frac{\partial L}{\partial \dot q_{i}}) - \frac{\partial L}{\partial q_{i}} = Q_i, \text{ } i = 1 \dots n $$

Sitä, miksi liike- ja potentiaalienergian erotuksen dervioinnilla monella eri tavoin, on jotain järkeä, ei voi ymmärtää ylläolevaa yhtälöä tuijottamalla eikä "maalaisjärjellä." Yritin ymmärtää yhtälön johtamista ja mielenkiintoiselta vaikutti, mutta syvimpiä perusteita en osaa selittää. Variaatiolaskenta oli keskeistä ja sitä en lähde tässä esittelemään.

Ulkoisten voimien laskemista selitetään mm wikipedia artikkelissa generalized forces. Löysin samalle asialle mielestäni yksinkertaisemman muotoilun. Sen mukaan

$$ Q_i = \frac{\partial \dot W}{\partial \dot q_i} \text{ ja } \dot W = \mathbf F\cdot \mathbf v$$

Voimien $Q_i$ yhtälö on johdettu virtuaalityön käsitteen avulla. Virtuaalityön idea on laskea järjestelmän siirtämiseen tilasta A tilaan B tarvittava energia kahdella toisistaan äärettömän vähän poikkeavaa polkua pitkin. Jonkun teoreeman mukaan virtuaalityö minimoituu fysiikan lakeja vastaavalla polulla. Tai niin ymmärsin teorioita äkkiseltään silmäilemällä. Toisaalta $\mathbf F\cdot \mathbf v$ on tehon kaava. Lagrangen mekaniikka on tunnettu siitä, että yhtälöiden fysikaalinen tulkinta "maalaisjärjellä" ei ola aina helppoa. Jätetään siis enempi hämmästely.

Yhtälö näyttää jokatapauksessa pahalta. Onneksi voimme jättää sen ratkaisemisen tietokoneelle, kunhan ensin onnistumme kirjoittamaan nosturin liike- ja potentiaalienergian lausekkeet. Siinäkin suurimman osan työstä jätämme tietokoneelle.

Koordinaatisto

Lagrangen mekaniikan mukaisissa esimerkkiratkaisuissa neuvottiin valitsemaan koordinaatisto niin, että järjestelmän liike- ja potentiaalienergian lausekkeiden kirjoittaminen on mahdollisimman helppoa. Eli luovutaan xy-koordinaatistosta ja kuvataan vaunun ja kuorman paikka joidenkin muiden muuttujien avulla.

Tämä wikipedia artikkeli auttanee ymmärtämään koordinaatiston valintaa. (Ei se kyllä helppoa ollut.)

Monessa webistä löytämässäni esimerkissä oli heiluri ja sen tilaa kuvaaviksi muuttujiksi oli yleensä valittu heilurin varren pituus ja sen poikkeama pystysuunnasta. Tein samoin. Valitsin koordinaatit niin, että $x$ on massan $m_{car}$ liikkeen suuntainen koordinaatti, koordinaatti $\theta$ kertoo vaijerin poikkeaman pystysuunnasta radiaaneina ja $L$ on vaijerin pituuden suuntainen koordinaatti. Näin voimat vaikuttavat koordinaattien suuntaisesti, mikä tekee yhtälöiden kirjoittamisen helpoksi.

Jatkossa käytän "tavallisesta" koordinaatistosta nimeä "xy-koordinaatisto" ja Lagrangen mekaniikan vaatimasta koordinaatistosta nimeä "q-koordinaatisto", koska valitsin (yleisen käytännön mukaisesti) q:n vektoriksi, joka kertoo järjestelmän osien paikan valitsemassani koordinaatistossa.

$$q(t) = \begin{bmatrix} x_{car}(t) \\ \theta(t) \\ L(t) \end{bmatrix}$$

Olenko ymmärtänyt oikein?

Sovellan Lagrangen mekaniikkaa ensimmäistä kertaa, joten vahvin todiste siitä, että olen ymmärtänyt asiat oikein ovat järkevän näköiset tulokset, joita lopulta sain aikaiseksi. Matemaattisesti ottaen tuo ei vielä todista mitään, mutta minä olen tyytyväinen.

Lagrangen yhtälön soveltaminen on insinöörilaskentoa, mutta sen johtaminen ja johtamiseen perehtyminen matematiikkaa. Tällä kertaa tyydyin verryttelemään insinöörilaskentoa. Jos olisin ratkaisemassa jotain lentokoneen tai laivan ohjaamiseen liittyvää ongelmaa, kävisin juttelemassa jonkun matemaatikon kanssa.

Koodin tuottaminen simulointiin ja optimointiin

Aletaan laatia python-kielistä ohjelmaa, joka kirjoittaa python-kieliset funktiot, joita tarvitaan siltanosturin simulointiin ja optimaalisen ohjauksen laskemiseen.

Tilayhtälöiden johtaminen

Järjestemän paikan q-koordinaatistossa ilmoittavat vaunun paikka, nostovaijerin pituus ja kulma, jonka se muodostaan pystysuoran kanssa.

Järjestelmän tilamuuttujiksi tarvitaan myös nopeudet, eli paikkavektorin derivaatta ajan suhteen.

Kun tilayhtälöitä johdetaan Lagrangen mekaniikan mukaisesti, tarvitaan vielä vektorin q toinen derivaatta ajan suhteen.

In [21]:
x_car, theta, L, t = sp.symbols('x_car, theta, L, t')
q = sp.Matrix([x_car(t), theta(t), L(t)])
sysdim = len(q)
display("paikka ",q)

dqdt = sp.Matrix(sp.diff(q,t))
display("nopeus ",dqdt)

d2qdt2 = sp.Matrix(sp.diff(q,t,2))
display("kiihtyvyys ",d2qdt2)
'paikka '
$$\left[\begin{matrix}\operatorname{x_{car}}{\left (t \right )}\\\theta{\left (t \right )}\\L{\left (t \right )}\end{matrix}\right]$$
'nopeus '
$$\left[\begin{matrix}\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\\\frac{d}{d t} \theta{\left (t \right )}\\\frac{d}{d t} L{\left (t \right )}\end{matrix}\right]$$
'kiihtyvyys '
$$\left[\begin{matrix}\frac{d^{2}}{d t^{2}} \operatorname{x_{car}}{\left (t \right )}\\\frac{d^{2}}{d t^{2}} \theta{\left (t \right )}\\\frac{d^{2}}{d t^{2}} L{\left (t \right )}\end{matrix}\right]$$

Järjestelmän potentiaali- ja liike-energia

Lagrangen mekaniikassa tilayhtälöiden johtamisessa käytetään apuna järjestelmän potentiaali- ja liike-energiaa. Ellei järjestelmään vaikuttaisi ulkoisia ohjaus- tai kitkavoimia, voimia ei tarvitsi miettiä lainkaan.

järjestelmän potentiaali- ja liike-energia on helppo esittää xy-koordinaatistossa, joten laskemme vaunun ja kuorman paikat ja nopeudet xy-koordinaatistossa q-koordinaattien avulla. Vaunun osalta ne tosin ovat samat.

Kuorman paikka xy-koordinaatistossa

(Tässä kohtaa sopii kerrata tai esitellä trigonometriaa.)

Tehdään vektoreista ja matriiseista sympy-modulin toivomia käyttämällä sp.Matrix-funktiota.

In [22]:
xy_load = sp.Matrix([x_car(t) + L(t)*sp.sin(theta(t)), -L(t)*sp.cos(theta(t))])
xy_load
Out[22]:
$$\left[\begin{matrix}L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} + \operatorname{x_{car}}{\left (t \right )}\\- L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )}\end{matrix}\right]$$

Kuorman potentiaalienergia.

Potentiaalienergia on kappaleen massa kerrottuna maan vetovoiman kiihtyvyydellä ja korkeudella vertailupisteeksi valitusta 0-tasosta: $W = mgh$ 0-tasoksi on valittu kuvassa 1 ylhäällä vaunun kiskon tasossa, joten tässä tapauksessa potentiaalienergia on aina negatiivinen. Vähän hämäävää, mutta teki yhtälöistä yksinkertaisia.

Tuo tieto riittää, mutta minua jäi vaivaamaan mistä potentiaalienergian kaava tulee?

Vapaasti putoavan kappaleen kiihtyvyys Maan vetovoimakentässä on $g=9.81$. Koska $F = ma$ putoavaan kappaleeseen vaikuttaa siis voima $F_{g} = mg$. Kappale pysyy paikallaan jos siihen kohdistuu maan vetovoiman kumoava vastakkaissuuntainen voima $F_{tuki} = mg$, jota yleensä kutsumme kappaleen painoksi.

Energia on yhtä kuin tehty työ. Potentiaalienergia on kappaleen nostamiseen tehty työ tai sen putoamisessa vapautuva energia. Työ on $W = F\cdot \Delta ds$ eli kappaleen nostamiseen tehty työ on $ W = mg\cdot \Delta s$

Jos voima riippuu paikasta, kertolasku ei riitä vaan työ lasketaan integroimalla: $W = \int F(s)\cdot ds$

Arvion — approksimaation — tehdylle työlle saisi laskemalla yhteen pienen pienillä matkoilla tehdyt työt: $W = F(0)ds + F(ds)ds + F(2*ds)ds \dots$ Mitä lyhyemmäksi ds:n valitsee, sitä tarkempi approksimaatiosta tulee. Kun ds lähestyy ääretöntä, summasta tulee integraali.

Potentiaalienergian lauseke pythonilla

Pythonissa indeksointi aloitetaan 0:sta, joten xy_load[1] tarkoittaa xy_load:in toista alkiota eli y-komponenttia eli kuorman korkeutta.

In [23]:
M_load, g = sp.symbols('M_load, g')
Vg = M_load * g * xy_load[1]
display('kuorman korkeus ',xy_load[1],'potentiaalienergia ',Vg)
'kuorman korkeus '
$$- L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )}$$
'potentiaalienergia '
$$- M_{load} g L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )}$$

Liike-energia

Liike-energia on kappaleen massa kerrottuna nopeuden neliöllä ja jaettuna kahdella: $W = \frac{1}{2}m\cdot v^2$

Tuon kaavan tietäminen riittää, mutta minua alkoi taas vaivata, mistä kaava tulee. Millä perusteella se on tuon muotoinen?

Jos kappaleen potentiaalienergia ei muutu, voiman tekemä työ kuluu kappaleen nopeuden kiihdyttämiseen eli muuttuu liike-energiaksi. Siirretään kappaletta paikasta 0 paikkaan $s_f$ voimalla $F(s)$ Ei ole väliä, millaista rataa kappaletta siirretään, joten tyydytään suoraviivaiseen liikkeeseen.

$W = \int_0^{s_f} F(s)ds$ koska työ on voima kertaa matka.

$W = \int_0^T F(s(t))v(t)\,dt = m\int_0^T \frac{d}{dt}v(t)v(t)\,dt$, koska $ds = vdt, F = ma$ ja $a(s(t))$ voidaan korvata $a(t)$:llä, koska voiman tekemä työ on sama siitä riippumatta, mitä rataa kuljetaan (myönnetään, tämä ei selitä asiaa kunnolla.)

Integraali näyttää siltä, että se ratkeaisi suoraan jollain integrointisäännöllä. Niitä on kuitenkin valtavasti, en muista ensimmäistäkään eikä minulla ole enää tallella kaavakokoelmaa, josta ratkaisu löytyisi.

Katsotaan, mihin muotoon sympy kirjoittaa W:n lausekkeen:

In [24]:
v, s, t, T = sp.symbols('v, s, t, T')
W = m * sp.integrate(sp.diff(v(t), t) * v(t), (t, 0, T))
W.subs(v(0), 0)  # alkunopeus = 0
Out[24]:
$$\frac{m}{2} v^{2}{\left (T \right )}$$

Liike-energia xy-koordinaatistossa

xy-koordinaatistossa voimme esittää nopeuden komponenttimuodossa

In [25]:
v_x, v_y = sp.symbols('v_x, v_y')
v = sp.Matrix([v_x, v_y])
v
Out[25]:
$$\left[\begin{matrix}v_{x}\\v_{y}\end{matrix}\right]$$

Pythagoraan teoreeman mukaan $v² = v_x^2 + v_y²$, joten Pythagorasta voisi muistella tässä kohdin. Nopeuden neliön saamme lsakemalla nopeusvektorin pistetulon itsensä kanssa.

In [26]:
m = sp.symbols('m')
K = sp.simplify(1 / 2 * m * v.dot(v))
K
Out[26]:
$$0.5 m \left(v_{x}^{2} + v_{y}^{2}\right)$$

Vaunun nopeus ja liike-energia

In [27]:
M_car = sp.symbols('M_car')
v_car = sp.diff(x_car(t), t)
T1 = (M_car * v_car**2) / 2
display(v_car, T1)
$$\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}$$
$$\frac{M_{car}}{2} \left(\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right)^{2}$$

Kuorman nopeus xy-koordinaatistossa

Kuorman liike-energia on helppo esittää xy-koordinaatistossa, joten tarvitsemme kuorman nopeuden xy-koordinaatistossa esitettynä q-koordinaatiston suureiden avulla. Kulma $\theta$ ja vaijerin pituus kertovat, paljonko kuorma on edellä tai jäljessä nosturin vaunua x-suunnassa ja paljonko kuorma on nosturin kiskon alapuolella y-suunnassa.

In [28]:
xy_load
Out[28]:
$$\left[\begin{matrix}L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} + \operatorname{x_{car}}{\left (t \right )}\\- L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )}\end{matrix}\right]$$

Kuorman nopeus on sen paikan derivaatta ajan suhteen. Derivoimme koko matriisin kerrallaan.

In [29]:
v_load = sp.Matrix(sp.diff(xy_load, t))
v_load
Out[29]:
$$\left[\begin{matrix}L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} + \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} + \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\\L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} - \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\end{matrix}\right]$$

Kuorman liike-energia

Liike-energia = $\frac{1}{2}M_{load}\underline v \cdot \underline v $

In [30]:
v2_krm = v_load.dot(v_load)
v_krm = sp.sqrt(v2_krm)
T2 = sp.simplify(1 / 2 * M_load * v2_krm)
T2
Out[30]:
$$0.5 M_{load} \left(L^{2}{\left (t \right )} \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} + 2 L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + 2 \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + \left(\frac{d}{d t} L{\left (t \right )}\right)^{2} + \left(\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right)^{2}\right)$$

Jousen potentiaalienergia

Jousi?

Jousta ei nosturissa vielä ole, mutta se olisi helppo lisätä liikeyhtälöiden johtamiseen.

(Kirjoitin helppo, mutta se ei tarkoita, etteikö jousen lisääminen teettäisi työtä ja etteikö se saattaisi tuoda yllätyksiä mukanaan. Muutama ihan helppo idea on pitkän puurtamisen jälkeen osoittautunut huonoksi ideaksi. Usein "hyvän" idean typeryyden ymmärtää vasta yritettyään toteuttaa sitä.)

Nostovaijeri venyy, joten sen jatkona voisi ajatella olevan jousen. Niin kauan, kuin voima $F_{load}$ nostaa suoraan kuormaa, jousen lisäämisessä ei ole mieltä (kokeilin monta kertaa ennen kuin tajusin), mutta jos esimerkiksi mallintaisimme vaijeria kelaavan laitteiston hitauden, vaijerin jatkona oleva jousi toisi nosturin dynamiikkaan mielenkiintoisen lisäefektin. Lisäksi pääsisimme harjoittelemaan pyörivän liikkeen dynamiikan mallintamista.

Uusia tilamuuttujia olisivat kelan pyörähtämät kierrokset, kelan pyörimisnopeus ja jousen venymä $\Delta s$. Uusi parametri olisi kelan ja sitä pyörittävän laitteiston inertiamomemtti, joka vastaa suoraviivaisen liikkeen massaa. Voiman $F_{load}$ sijasta toisena ohjaussuurena olisi kelaa pyörittävä momentti. Vaijerin pituus laskettaisiin kelan kierroksista ja kuorman paikka olisi vaijerin pituus lisättynä $\Delta s$:llä.

Jousen puristumaan tai venymään $\Delta s$ varastoitunut potentiaalienergia on $W = \frac{1}{2} k \cdot \Delta s^{2}$, missä $k$ on jousen jäykkyyttä kuvaava jousivakio. Tämä potentiaalienergia lisättäisiin muihin potentiaalienergioihin seuraavassa vaiheessa.

Tilayhtälöiden johtaminen Lagrangen periaatteella

Lagrangen mekaniikkaa sovellettaessa muodostetaan langrangen funktio L, joka on järjestelmän kineettisen ja potentiaalienergian erotus ja ratkaistaan yhtälöistä

$$ \frac{d}{dt}(\frac{\partial L}{\partial \dot q_{i}}) - \frac{\partial L}{\partial q_{i}} = Q_i, \text{ } i= 1 \dots n $$ $$ Q_i = \frac{\partial \dot W}{\partial \dot q_i} \text{ ja } \dot W = \mathbf F\cdot \mathbf v$$

vaunun ja kuorman kiihtyvyydet $\ddot q(t) = $ {{d2qdt2}} niin, että saadaan yhtälö

$$\ddot q(t) = f(\dot q(t), q(t)) + F_u $$

Meidän siis pitää kirjoittaa järjestelmän potentiaali- ja liike-energian lausekkeet sekä lauseke ohjaus- ja kitkavoimien teholle.

Järjestelmään vaikuttavat ulkoiset voimat

Voimat, joilla järjestelmää ohjataan

Ohjausvoimat ovat vaunua siirtämään pyrkivä voima ja vaijeria kiristävä voima. $F_{car}$ vaikuttaa vaunun liikeradan suuntaan ja $F_{load}$ nostovaijerin suuntaan. Positiivinen $F_{car}$ vetää vaunua positiiviseen suuntaan. Positiivinen $F_{load}$ pyrkii lyhentämään vaijeria, siksi - merkki. (Tästä miinusmerkistä oli harmia monta kertaa.)

In [31]:
udim = 2  # apumuuttuja koodin generointia varten, ohjausmuuttujien määrä
F_car, F_load = sp.symbols('F_car, F_load')
W_u = F_car*v_car - F_load*sp.diff(L(t),t)
Q_u = [sp.diff(W_u,dq) for dq in dqdt]
Q_u
Out[31]:
$$\left [ F_{car}, \quad 0, \quad - F_{load}\right ]$$

Koska mikään ulkoinen voima tai momentti ei väännä vaijerin kulmaa, keskimmäinen ohjausvoima on nolla. Tulos on muutenkin itsestäänselvän näköinen. Hyvä niin.

Kitkavoimat

Oletin vaunun ja kuorman kitkat ja ilmanvastukset suoraan verrannollisiksi nopeuteen $F_{fr}= −Cf⋅v_F$, jotta sain ne sovitettua Lagrangen mekaniikan mukaiseen ratkaisuun. Pienillä nopeuksilla se ei ole ihan paha oletus. (Sain käsityksen, että nopeuden neliöön verrannollista kitkavoimaa ei voi ujuttaa Langrange mekaniikkaan.)

Kitkavoima jarruttaa matkan tekoa, joten se on erimerkkinen kuin nopeus.

Nostovaijeriin kelaamiseen liittyvän kitkan suuntaa on syytä hetki miettiä. Kun nostovaijerin nopeus on positiivinen, vaijeri pitenee ja pitenemistä jarruttava voima osoittaa vaunun suuntaan, jonka olen valinnut positiiviseksi nostavan voiman suunnaksi. Kitkan tekemän virtuaalisen tehon etumerkiksi lienee siis laitettava plus.

In [32]:
cfr_car, cfr_load, cfr_kela = sp.symbols('cfr_car, cfr_load, cfr_kela')
W_fr = -cfr_car * v_car * v_car - cfr_load * v_load.dot(
    v_load) + cfr_kela * sp.diff(L(t), t) * sp.diff(L(t), t)
Q_fr = [sp.diff(W_fr, dq) for dq in dqdt]
Q_fr
Out[32]:
$$\left [ - 2 cfr_{car} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} - cfr_{load} \left(2 L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} + 2 \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} + 2 \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right), \quad - cfr_{load} \left(2 \left(L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} - \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\right) L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} + 2 \left(L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} + \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} + \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right) L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )}\right), \quad 2 cfr_{kela} \frac{d}{d t} L{\left (t \right )} - cfr_{load} \left(- 2 \left(L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} - \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\right) \cos{\left (\theta{\left (t \right )} \right )} + 2 \left(L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} + \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} + \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right) \sin{\left (\theta{\left (t \right )} \right )}\right)\right ]$$

Lagrangen mekaniikan mukaan järjestelmän liike- ja potentiaalienergian erotukselle $L$ ja ulkoisille voimille $\mathbf Q$ pätee:

$$ \frac{d}{dt}(\frac{\partial L}{\partial \dot q_{i}}) - \frac{\partial L}{\partial q_{i}} = Q_i, \text{ } i= 1 \dots n $$

Pahalta näyttää, mutta onneksi derivaattojen lausekkeet voi kirjoittaa suoraan sp.diff-kommennoiksi. Annetaan ohjelman laskeskella ja katsotaan mitä tulee.

Seuraavassa käytän L:n sijasta muuttujaa LE, koska L on jo varattu nostovaijerin pituudeksi.

In [33]:
LE = T1 + T2 - Vg 
LE
Out[33]:
$$\frac{M_{car}}{2} \left(\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right)^{2} + M_{load} g L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} + 0.5 M_{load} \left(L^{2}{\left (t \right )} \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} + 2 L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + 2 \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + \left(\frac{d}{d t} L{\left (t \right )}\right)^{2} + \left(\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right)^{2}\right)$$
In [34]:
display(q,dqdt, d2qdt2)
$$\left[\begin{matrix}\operatorname{x_{car}}{\left (t \right )}\\\theta{\left (t \right )}\\L{\left (t \right )}\end{matrix}\right]$$
$$\left[\begin{matrix}\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\\\frac{d}{d t} \theta{\left (t \right )}\\\frac{d}{d t} L{\left (t \right )}\end{matrix}\right]$$
$$\left[\begin{matrix}\frac{d^{2}}{d t^{2}} \operatorname{x_{car}}{\left (t \right )}\\\frac{d^{2}}{d t^{2}} \theta{\left (t \right )}\\\frac{d^{2}}{d t^{2}} L{\left (t \right )}\end{matrix}\right]$$
In [35]:
#dl_1 = [sp.diff(LE, dqdt[i]) for i in range(0, sysdim)]
dl_1 = [sp.diff(LE, dq) for dq in dqdt]
dl_1
Out[35]:
$$\left [ M_{car} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + 0.5 M_{load} \left(2 L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} + 2 \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} + 2 \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right), \quad 0.5 M_{load} \left(2 L^{2}{\left (t \right )} \frac{d}{d t} \theta{\left (t \right )} + 2 L{\left (t \right )} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right), \quad 0.5 M_{load} \left(2 \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + 2 \frac{d}{d t} L{\left (t \right )}\right)\right ]$$
In [36]:
dl_2 = [sp.diff(LE, qi) for qi in q]
dl_2
Out[36]:
$$\left [ 0, \quad - M_{load} g L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} + 0.5 M_{load} \left(- 2 L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + 2 \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right), \quad M_{load} g \cos{\left (\theta{\left (t \right )} \right )} + 0.5 M_{load} \left(2 L{\left (t \right )} \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} + 2 \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \theta{\left (t \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\right)\right ]$$
In [37]:
syseqs = [sp.diff(dl_1[i], t) - dl_2[i] - (Q_u[i]+Q_fr[i]) for i in range(0, sysdim)]
In [38]:
kiihtyvyydet = sp.simplify(sp.solve(syseqs, d2qdt2))
kiihtyvyydet
Out[38]:
$$\left \{ \frac{d^{2}}{d t^{2}} L{\left (t \right )} : - \frac{F_{car}}{M_{car}} \sin{\left (\theta{\left (t \right )} \right )} - \frac{F_{load}}{M_{load}} - \frac{F_{load}}{M_{car}} \sin^{2}{\left (\theta{\left (t \right )} \right )} + g \cos{\left (\theta{\left (t \right )} \right )} + L{\left (t \right )} \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} + \frac{2.0 cfr_{kela}}{M_{load}} \frac{d}{d t} L{\left (t \right )} - \frac{2.0 cfr_{load}}{M_{load}} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} - \frac{2.0 cfr_{load}}{M_{load}} \frac{d}{d t} L{\left (t \right )} + \frac{2.0 cfr_{car}}{M_{car}} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + \frac{2.0 cfr_{kela}}{M_{car}} \sin^{2}{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}, \quad \frac{d^{2}}{d t^{2}} \theta{\left (t \right )} : \frac{1}{M_{car} M_{load} L{\left (t \right )}} \left(- F_{car} M_{load} \cos{\left (\theta{\left (t \right )} \right )} - 0.5 F_{load} M_{load} \sin{\left (2.0 \theta{\left (t \right )} \right )} - M_{car} M_{load} g \sin{\left (\theta{\left (t \right )} \right )} - 2.0 M_{car} M_{load} \frac{d}{d t} L{\left (t \right )} \frac{d}{d t} \theta{\left (t \right )} - 2.0 M_{car} cfr_{load} L{\left (t \right )} \frac{d}{d t} \theta{\left (t \right )} - 2.0 M_{car} cfr_{load} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + 2.0 M_{load} cfr_{car} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + M_{load} cfr_{kela} \sin{\left (2.0 \theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\right), \quad \frac{d^{2}}{d t^{2}} \operatorname{x_{car}}{\left (t \right )} : \frac{1}{M_{car}} \left(F_{car} + F_{load} \sin{\left (\theta{\left (t \right )} \right )} - 2.0 cfr_{car} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} - 2.0 cfr_{kela} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\right)\right \}$$

Siinä meillä on tilayhtälö toisen kertaluvun differentiaaliyhtälöinä. Numeerinen ratkaisu vaatii tilayhtölän esittämistä ensimmäisen kertaluvun differentiaaliyhtlöinä, joten alamme vääntää ylläolevaa tilayhtälöä muotoon $\dot{\underline{z}}(t) = \underline{f}(\underline{z}(t))$

In [39]:
f = sp.Matrix([0 for i in range(0, 2 * sysdim)])
for i in range(0, sysdim):
    f[i] = kiihtyvyydet[d2qdt2[i]]
    f[i + 3] = dqdt[i]
display(f)
$$\left[\begin{matrix}\frac{1}{M_{car}} \left(F_{car} + F_{load} \sin{\left (\theta{\left (t \right )} \right )} - 2.0 cfr_{car} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} - 2.0 cfr_{kela} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\right)\\\frac{1}{M_{car} M_{load} L{\left (t \right )}} \left(- F_{car} M_{load} \cos{\left (\theta{\left (t \right )} \right )} - 0.5 F_{load} M_{load} \sin{\left (2.0 \theta{\left (t \right )} \right )} - M_{car} M_{load} g \sin{\left (\theta{\left (t \right )} \right )} - 2.0 M_{car} M_{load} \frac{d}{d t} L{\left (t \right )} \frac{d}{d t} \theta{\left (t \right )} - 2.0 M_{car} cfr_{load} L{\left (t \right )} \frac{d}{d t} \theta{\left (t \right )} - 2.0 M_{car} cfr_{load} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + 2.0 M_{load} cfr_{car} \cos{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + M_{load} cfr_{kela} \sin{\left (2.0 \theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\right)\\- \frac{F_{car}}{M_{car}} \sin{\left (\theta{\left (t \right )} \right )} - \frac{F_{load}}{M_{load}} - \frac{F_{load}}{M_{car}} \sin^{2}{\left (\theta{\left (t \right )} \right )} + g \cos{\left (\theta{\left (t \right )} \right )} + L{\left (t \right )} \left(\frac{d}{d t} \theta{\left (t \right )}\right)^{2} + \frac{2.0 cfr_{kela}}{M_{load}} \frac{d}{d t} L{\left (t \right )} - \frac{2.0 cfr_{load}}{M_{load}} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} - \frac{2.0 cfr_{load}}{M_{load}} \frac{d}{d t} L{\left (t \right )} + \frac{2.0 cfr_{car}}{M_{car}} \sin{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )} + \frac{2.0 cfr_{kela}}{M_{car}} \sin^{2}{\left (\theta{\left (t \right )} \right )} \frac{d}{d t} L{\left (t \right )}\\\frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}\\\frac{d}{d t} \theta{\left (t \right )}\\\frac{d}{d t} L{\left (t \right )}\end{matrix}\right]$$

Joitakuita saattoi yllä hämmentää lauseke f[i] = kiihtyvyydet[d2qdt2[i]].

sp.solve antaa yhtälön ratkaisun niin kutsuttuna dictionaryna, josta alla r on esimerkki. Dictionaryn alkiohin ei voi viitata numeroilla vaan pitää käyttää "avaimia", joita r:ssä ovat 'a' ja 'b'

In [40]:
r = {'a': 5, 'b': 7}
eka = r['a']
eka
Out[40]:
$$5$$

Yhtälöryhmän ratkaisussa avaimina ovat muuttujien toiset derivaatat, joten voin käyttää avaimina matriisin d2qdt2 alkioita {{d2qdt2}} samoin kuin käytän alla listan ab alkioita.

(Itse ainakin hämmästelen, miten keksinkään tällaisen ratkaisun ;-)

In [41]:
ab = ['a', 'b']
eka = r[ab[0]]
eka
Out[41]:
$$5$$

Koodin generointi simulointia varten

Alla ensiksi pari esimerkkiä selittämään, mitä kannattaa tehdä ennen koodin generoimista. Emme esimerkiksi halua python-koodiin turhaa Matrix-määrittelyä. (Tämä on amatöörin ajattelua, kokeneempi koodari tietänee paremman ratkaisun.)

In [42]:
esim1 = sp.Matrix([1, 2, 3])
print("esim1 tuottaisi tällaista koodia: ", esim1)
esim1
esim1 tuottaisi tällaista koodia:  Matrix([[1], [2], [3]])
Out[42]:
$$\left[\begin{matrix}1\\2\\3\end{matrix}\right]$$
In [43]:
esim2 = [x for x in esim1]
print("esim2 tuottaisi tällaista koodia: ", esim2)
esim2
esim2 tuottaisi tällaista koodia:  [1, 2, 3]
Out[43]:
$$\left [ 1, \quad 2, \quad 3\right ]$$

Systeemissä on vakioparametrejä, joita on pikkuisen hankala käsitellä koodin generoinnissa. Luon omaan moduulinsa olion par, jonka elementteinä ovat järjestelmän parametrit.

Syystä, joka selviää myöhemmin, upotan aluksi vakioiden nimet lausekkeen pp() sisään.

In [44]:
par_subs = []

params = [
    ('sysdim', 3), ('udim', 2), ('g', 9.81), ('M_car', 1.0), ('M_load', 3.0),
    ('T_f', 6.0), ('x0_car', 0.0), ('xf_car', 12.0), ('L_0', 8.0),
    ('L_f', 4.0), ('cfr_car', 0.1), ('cfr_load', 0.5), ('cfr_kela', 0.0),
    ('C_W', 1.0), ('C_F1', 10.0), ('C_F2', 2.0), ('u1lim', 100.0),
    ('u2lim', 100.0)
    # ('cf_1', 1000.0),
    # ('cf_2', 1000.0),
    # ('cf_3', 1000.0),
    # ('cf_4', 1000.0),
    # ('cf_5', 10000.0),
    # ('cf_6', 1000.0)    
]

par_subs = [(Name, "pp(" + str(Name) + ")")
            for i, (Name, Value) in enumerate(params)]
par_subs
Out[44]:
[('sysdim', 'pp(sysdim)'),
 ('udim', 'pp(udim)'),
 ('g', 'pp(g)'),
 ('M_car', 'pp(M_car)'),
 ('M_load', 'pp(M_load)'),
 ('T_f', 'pp(T_f)'),
 ('x0_car', 'pp(x0_car)'),
 ('xf_car', 'pp(xf_car)'),
 ('L_0', 'pp(L_0)'),
 ('L_f', 'pp(L_f)'),
 ('cfr_car', 'pp(cfr_car)'),
 ('cfr_load', 'pp(cfr_load)'),
 ('cfr_kela', 'pp(cfr_kela)'),
 ('C_W', 'pp(C_W)'),
 ('C_F1', 'pp(C_F1)'),
 ('C_F2', 'pp(C_F2)'),
 ('u1lim', 'pp(u1lim)'),
 ('u2lim', 'pp(u2lim)')]

Valmistaudutaan generoimaan simuloinnissa tarvittava järjestelmän tilan derivaatan palauttava funktio. Muuttujien vaihdoilla saamme kuusi ensimmäisen kertaluvun differentiaalyhtälöä. Esimerkiksi:

  • vaunun paikan toinen derivaatta ajan suhteen --> vaunun nopeuden derivaatta: $ \frac{d^{2}}{dt^{2}}x_{car}(t)\mapsto \frac{d}{dt}v_{car}(t)$
  • vaunun paikan derivaatta ajan suhteen --> vaunun nopeus: $\frac{d}{dt}x_{car}(t) \mapsto v_{vnu}$

Tästä eteenpäin emme tarvitse tietoa muuttujien aikariippuvuudesta, joten hävitämme niistä aika-argumentin sijoituksilla $L(t) \mapsto L, v_{car}(t) \mapsto v_{vnu},\dots$

Tein korvauksen $v_{car}(t) \mapsto v_{vnu}$ koska muuttuja $v_{car}(t)$ on jo käytössä ja saattaa vielä esiintyä jossain lausekkeessa. (Tämä on pikkuisen epämääräinen temppu. Ehkä joskus korjaan tämän.)

In [45]:
omega, v_L, v_vnu = sp.symbols('omega v_L v_vnu')
dx_subs = [(sp.diff(x_car(t), t), v_vnu), (sp.diff(theta(t), t), omega),
           (sp.diff(L(t), t), v_L), (theta(t), theta), (x_car(t), x_car),
           (L(t), L)]
dx_subs
Out[45]:
$$\left [ \left ( \frac{d}{d t} \operatorname{x_{car}}{\left (t \right )}, \quad v_{vnu}\right ), \quad \left ( \frac{d}{d t} \theta{\left (t \right )}, \quad \omega\right ), \quad \left ( \frac{d}{d t} L{\left (t \right )}, \quad v_{L}\right ), \quad \left ( \theta{\left (t \right )}, \quad \theta\right ), \quad \left ( \operatorname{x_{car}}{\left (t \right )}, \quad x_{car}\right ), \quad \left ( L{\left (t \right )}, \quad L\right )\right ]$$
In [46]:
dx0 = f.subs(dx_subs)
dx = [dxi for dxi in dx0]

Järjestelmän tilavektori ja ohjausvektori (jälkimmäiseen ympätty muutama parametri)

In [47]:
yy = [v_vnu, omega, v_L, x_car, theta, L]
yy
Out[47]:
$$\left [ v_{vnu}, \quad \omega, \quad v_{L}, \quad x_{car}, \quad \theta, \quad L\right ]$$
In [48]:
uu = [F_car, F_load, M_car, M_load, cfr_car, cfr_load, cfr_kela]

Simulaattorissa vaijerin voi "lukita". Lukittaessa vaijerin pituuden muutosnopeus asetetaan nollaksi. Lisäksi voima asetetaan kuorman painon suuruiseksi, vaikkei sillä ole kuin kosmeettista merkitystä. Tyylikkäämpää olisi laskea, mikä todellinen voima vaijeria kuormittaa lukitustilanteessa. Täytyy miettiä, onko sitä vaikea laskea.

In [49]:
F_lukittu = M_load*g

Numeerinen integrointi

Kirjoitetaan simuloinnissa tarvittava funktio fdxdt(), joka palauttaa järjestelmämuuttujien derivaatat järjestelmän tilan ja ohjausten funktiona. Videopelimäisessä simulaattorissa tätä funktiota käytetään seuraavasti: Generoimamme funktion fdxdt() nimi annetaan argumentiksi integrointifuntiolle rk4, joka kutsuu sitä laskiessaan järjestelmän siirtymää $yy(t) \to yy(t+\Delta t)$

from dxdt_sympy import fdxdt
from rk4 import rk4

def liike():
    yy = [vaunu.v, kuorma.w, kuorma.Lv, vaunu.x, kuorma.theta, kuorma.L]
    uu = [vaunu.F, kuorma.F, vaunu.M, kuorma.M, vaunu.Cf, kuorma.Cf]
    yy = rk4(fdxdt, yy, uu, Aika.Dt)
    [vaunu.v, kuorma.w, kuorma.Lv, vaunu.x, kuorma.theta, kuorma.L] = yy

rk4 taas on seuraavanlainen järjestelmää yhden aika-askeleen eteenpäin siirtävä funktio (runge-kutta integrointi)

def rk4(fdydt, yy, uu, dt):
    kk1 = fdydt(yy, uu)
    yk1 = [yy[i] + k1 * dt / 2.0 for i, k1 in enumerate(kk1)]
    kk2 = fdydt(yk1, uu)
    yk2 = [yy[i] + k2 * dt / 2.0 for i, k2 in enumerate(kk2)]
    kk3 = fdydt(yk2, uu)
    yk3 = [yy[i] + k3 * dt for i, k3 in enumerate(kk3)]
    kk4 = fdydt(yk3, uu)
    yy1 = [
        y + dt / 6.0 * (kk1[i] + 2.0 * kk2[i] + 2.0 * kk3[i] + kk4[i])
        for i, y in enumerate(yy)
    ]
    return yy1

Symbolisista lausekkeista merkkijonoja

Edelläjohdetut lausekkeet vaativat pientä manipulointia, ennenkuin ne voidaan kirjoittaa tiedostoon python-ohjelmaksi. Osa lausekkeista on listoja, jotka pitää käsitellä alkio kerrallaan niin, että tuloksenakin on lista.

Lausekkeesta expr syntyvä koodi kirjoitetaan tiedostoon fil, sitä sisennetään merkkijonon indent verran ja symbolisen lausekkeen osia korvataan listan par_subs mukaan.

In [50]:
def wr_expr(fil, expr, indent, par_subs):
    if (type(expr) == list):
        dim = len(expr) - 1
        fil.write('[')
        for i, ex in enumerate(expr):
            wr_expr(fil, ex, indent, par_subs)
            if (i < dim):
                fil.write(',')
        fil.write(']')

    else:
        subs_ex(fil, expr, indent, par_subs)

Alla subs_ex(), eli lausekkeen, joka ei ole lista, kirjoittaminen koodiksi.

Generoitavassa koodissa tilamuuttujat esitetään vektorina zz. Ennen tämän funktion kutsumista järjestelmän tilamuuttujat on sympy.subs-funktiolla koodattu muotoon zz(i), koska en sympy.subs-komennolla onnistunut saamaan niitä muotoon zz[i]. Esimerkiksi $v_{car} = zz(0)$, mutta haluamme sen muotoon $zz[0]$

Seuraavassa lauseke muutetaan ensin merkkijonoksi ja sen jälkeen käytetään regular expressions komentoja, jotka ovat muotoa re.sub(r'korvattava merkkijono', r'millä korvataan', merkkijono)

r'[0-9]' tarkoittaa mitä hyvänsä numeroa välitltä 0-9. r'[0-9]*' tarkoittaa edellinen yksi tai useampi kertaa, eli mikä hyvänsä numerojono.

Edellinen käärittynä kaarisulkeisiin — r'([0-9]*)' — tarkoittaa, että löytyvään merkkijonoon voidaan myöhemmin viitata merkinnällä '\1'.

Kenoviiva — esimerkiksi '\(' — tarkoittaa, että kenoviivaa seuraava merkki — edellä '(' — otetaan sellaisenaan.

In [51]:
tst1 = re.sub(r'zz\(([0-9]*)\)', r'zz[\1]', "pp(M_load)*zz(8)")
tst1
Out[51]:
'pp(M_load)*zz[8]'
In [52]:
tst2 = re.sub(r'pp\(([^)]*)\)', r'par.\1', tst1)
tst2
Out[52]:
'par.M_load*zz[8]'
In [53]:
def subs_ex(fil, ex, indent, par_subs):
    try:
        ex = ex.subs(par_subs)
    except:
        pass

#    s0 = re.sub(r'\^', r'',str(ex))
    s1 = re.sub(r'zz\(([0-9]*)\)', r'zz[\1]', str(ex))
    s2 = re.sub(r'pp\(([^)]*)\)', r'par.\1', s1)
    s = re.sub(r'cf\(([0-9]*)\)', r'par.cf_\1', s2)
    fil.write(
        fill(
            s,
            initial_indent=indent,
            subsequent_indent="    ",
            break_long_words=False))

Kirjoitetaan simulointiin tarvittava systeemiyhtälö

In [54]:
with open('dxdt.py', 'w') as fil:
    fil.write("# -*- coding: utf-8 -*-\n")
    fil.write("from math import cos, sin \n")
    fil.write("from opt_params import par \n\n")
    fil.write("g = 9.81 \n\n")

    fil.write("def fdxdt(yy, uu):\n    ")
    #    all_params()
    wr_expr(fil, yy, "    ", [])
    fil.write(" = yy\n    ")
    wr_expr(fil, uu, "", [])
    fil.write(" = uu\n")
    fil.write("    dxdt =")
    wr_expr(fil, dx, " ", [])
    fil.write("\n")
    fil.write("    return dxdt\n")

    fil.write("#####################\n")

    fil.write("def Flukko(yy, uu):\n    ")
    wr_expr(fil, uu, "", [])
    fil.write(" = uu\n")
    fil.write("    return ")
    wr_expr(fil, F_lukittu, "    ", [])
    fil.write("\n")
In [55]:
# with open('dxdt.py', 'r') as fin:
#    print(fin.read())

Optimiohjauksen laskeminen

Järjestelmän alku- ja lopputila

Järjestelmä halutaan levosta lepoon, eli kaikki nopeudet halutaan nolliksi sekä alkutilassa yy_0 että lopputilassa yy_f. Vaunu halutaan siirtää paikasta x0_car paikkaan xf_car ja kuorma halutaan nostaa tai laskea korkeudesta -L_0 korkeuteen -L_f. Miinusmerkit tarvitaan, koska L kertoo vaijerin pituuden, eli kuinka paljon 0-tasoksi valitun vaunun alapuolella kuorma riippuu.

In [56]:
x0_car, xf_car, L_0, L_f = sp.symbols('x0_car, xf_car, L_0, L_f')
yy_0 = [0.0, 0.0, 0.0, x0_car, 0.0, L_0]
yy_f = [0.0, 0.0, 0.0, xf_car, 0.0, L_f]

Optimaalinen siirtymä alusta loppuun

Siirretään kuorma alkutilasta lopputilaan optimaalisesti laskemalla sellaiset vaunua ja kuormaa liikuttavat voimat, että kustannusfunktion $gl$ integraali siirtymän alkuhetkestä loppuhetkeen minimoituu ohjauksien $\underline u(t)$ suhteen.

$$ J = G_f(\underline x(T_f)) + \int_0^{T_f} gl(\underline x(t),\underline u(t),t) \,\mathrm{d}t $$

In [57]:
C_F1, C_F2, C_W, T_f = sp.symbols('C_F1, C_F2, C_W, T_f')
F_0 = M_load * g

gl = C_F1 * F_car**2 + C_F2 * (F_load - F_0 * (1.0 - (L_f - L_0) / T_f)
                               )**2 + C_W * (F_car * v_vnu - F_load * v_L)
#gl = C_F1*F_car**2 + C_F2*(F_load - F_0)**2 + C_W*(F_car*v_vnu - F_load*v_L)

Kustannuksessa sakotetaan neliöllisesti vaunua vetävän voiman ja kuormaa kannattelevan voiman poikkeamaa lepoarvosta. Lisäksi sakotetaan käytetystä tehosta: $C_{W}(F_{car}v_{car} - F_{load}v_{L})$ (Kuormaa nostavan voiman merkki on erisuuntainen kuin vaijerin kelausnopeus.)

Millä perusteella tällainen kriteeri? Se on helppo laskea ja vie kuorman paikasta toiseen turhia tempoilematta.

Vaunun siirtelyn ja vaijerin kelailun vaatimien tehojen summan lisäsin kustannukseen, koska minusta oli jännittävä nähdä, miten sellainen vaikuttaa. Onhan energian säästö järkevääkin.

Tehojen lausekkeesta voi arvata, että ilman neliöllisiä kustannuksia optimaaliset ohjaukset ovat nopeuksista riippuen $\pm\infty$. Ilman edes pientä neliöllistä kustannukset ohjaukset pitää rajata järkevälle välille. Tällöin lopputuloksena on nk bang-bang ohjaus, jossa ohjaukset hyppelevät laidasta laitaan. Käyttämälläni algoritmilla tällainen tehtävä ratkeaa huonosti. Sopivilla kertoimien $C_{W}, C_{F1} \text{ ja } C_{F2}$ suhteella pääsemme kuitenkin tulokseen, jossa tehon käyttö painottuu selvästi.

Numeerisesti yhtä vaikea tehtävä on heittää kuorma mahdollisimman nopeasti paikasta toiseen.

Loppukustannus

Optimointialgoritmin voi vaatia viemään järjestelmä tarkalleen lopputilaan. Tehtävä ratkesi niin muotoiltuna hyvin.

Kokeilin kuitenkin tarkan lopputilan vaatimisen sijasta nelilöllistä sakkoa poikkeamalle halutusta lopputilasta, joten dokumentoin senkin ratkaisun tässä samalla.

Lopputulos oli jokseenkin sama, mutta joskus numeerinen algoritmi antoi ratkaisun silloinkin, kun tarkalleen lopputilaan ei voinut päästä esimerkiksi siksi, että olin asettanut ohjausvoimille liian tiukat rajat.

In [58]:
cf = sp.symbols('cf')
G_f = [cf(i + 1) * (yy[i] - yy_f[i])**2 for i in range(0, 2 * sysdim)]
G_f
Out[58]:
$$\left [ v_{vnu}^{2} \operatorname{cf}{\left (1 \right )}, \quad \omega^{2} \operatorname{cf}{\left (2 \right )}, \quad v_{L}^{2} \operatorname{cf}{\left (3 \right )}, \quad \left(x_{car} - xf_{car}\right)^{2} \operatorname{cf}{\left (4 \right )}, \quad \theta^{2} \operatorname{cf}{\left (5 \right )}, \quad \left(L - L_{f}\right)^{2} \operatorname{cf}{\left (6 \right )}\right ]$$
In [59]:
lamda_f = [sp.diff(G_f[i], yy[i]) for i in range(0, 2 * sysdim)]
lamda_f
Out[59]:
$$\left [ 2 v_{vnu} \operatorname{cf}{\left (1 \right )}, \quad 2 \omega \operatorname{cf}{\left (2 \right )}, \quad 2 v_{L} \operatorname{cf}{\left (3 \right )}, \quad \left(2 x_{car} - 2 xf_{car}\right) \operatorname{cf}{\left (4 \right )}, \quad 2 \theta \operatorname{cf}{\left (5 \right )}, \quad \left(2 L - 2 L_{f}\right) \operatorname{cf}{\left (6 \right )}\right ]$$

Sovelletaan Pontryaginin minimiperiaatetta, joka johtaa nk kahden pisteen reuna-arvotehtävään. Sellaisten ratkaisemiseen on kehitetty numeerisia algoritmeja, joita varten väännämme ongelman niiden vaatimaan standardimuotoon.

Tilamuuttujien lisäksi tarvitaan yhtä monta nk. liittotilamuuttujaa.

H on Pontryaginin periaatteen mukainen Hamiltonin funktio ja yy_lt järjestelmän liittotila

In [60]:
# liittotilamuuttujat
lamda1, lamda2, lamda3, lamda4, lamda5, lamda6 = sp.symbols(
    'lamda_1, lamda_2, lamda_3, lamda_4, lamda_5, lamda_6')
yy_lamda = sp.Matrix([lamda1, lamda2, lamda3, lamda4, lamda5, lamda6])

# Hamiltonin funktio
H = yy_lamda.dot(dx) + gl

# Ohjausvektori
u = sp.Matrix([F_car, F_load])
u_lst = [F_car, F_load]

Optimaalinen ohjaus ajan funktiona löytyy minimoilla Hamiltonin funktio ohjauksen suhteen. Yksinkertaisen kustannusfunktiomme ansiosta minimi löytyy derivaatan 0-kohdasta

In [61]:
dHdu = [sp.diff(H, ui) for ui in u]
uopt = sp.simplify(sp.solve(dHdu, u))
uopt
Out[61]:
$$\left \{ F_{car} : \frac{0.5}{C_{F1} L M_{car}} \left(- C_{W} L M_{car} v_{vnu} - L \lambda_{1} + L \lambda_{3} \sin{\left (\theta \right )} + \lambda_{2} \cos{\left (\theta \right )}\right), \quad F_{load} : \frac{L_{0} M_{load}}{T_{f}} g - \frac{L_{f} M_{load}}{T_{f}} g + M_{load} g + \frac{0.5 C_{W}}{C_{F2}} v_{L} + \frac{0.5 \lambda_{3}}{C_{F2} M_{load}} - \frac{0.5 \lambda_{1}}{C_{F2} M_{car}} \sin{\left (\theta \right )} + \frac{0.5 \lambda_{3}}{C_{F2} M_{car}} \sin^{2}{\left (\theta \right )} + \frac{0.25 \lambda_{2}}{C_{F2} L M_{car}} \sin{\left (2.0 \theta \right )}\right \}$$

Hessian matriisia voidaan joissain tapauksissa käyttää apuna numeerisessa ratkaisussa.

HessH = [[0 for i in range(0,udim)]for i in range(0,udim)]
for i in range(0,udim):
    for j in range(0,udim):
        HessH[i][j] = sp.diff(dHdu[i],u[j])

Liittotilan tilayhtälöt (= liittotilamuuttujien aikaderivaatat) saadaan derivoimalla hamiltonin funktio varsinaisilla tilamuuttujilla.

In [62]:
dlt = sp.Matrix([-sp.simplify(sp.diff(H, yi)) for yi in yy])
dlt
Out[62]:
$$\left[\begin{matrix}- \frac{1}{L M_{car} M_{load}} \left(L M_{car} M_{load} \left(C_{W} F_{car} + \lambda_{4}\right) - 2.0 L M_{load} cfr_{car} \lambda_{1} - 2.0 L \lambda_{3} \left(M_{car} cfr_{load} - M_{load} cfr_{car}\right) \sin{\left (\theta \right )} - 2.0 \lambda_{2} \left(M_{car} cfr_{load} - M_{load} cfr_{car}\right) \cos{\left (\theta \right )}\right)\\- 2 L \lambda_{3} \omega - \lambda_{5} + \frac{2.0 cfr_{load}}{M_{load}} \lambda_{2} + \frac{2.0 \lambda_{2}}{L} v_{L}\\C_{W} F_{load} - \lambda_{6} - \frac{2.0 cfr_{kela}}{M_{load}} \lambda_{3} + \frac{2.0 cfr_{load}}{M_{load}} \lambda_{3} + \frac{2.0 cfr_{kela}}{M_{car}} \lambda_{1} \sin{\left (\theta \right )} - \frac{2.0 cfr_{kela}}{M_{car}} \lambda_{3} \sin^{2}{\left (\theta \right )} + \frac{2.0 \lambda_{2}}{L} \omega - \frac{cfr_{kela} \lambda_{2}}{L M_{car}} \sin{\left (2.0 \theta \right )}\\0\\- \frac{1}{L M_{car} M_{load}} \left(L M_{load} \lambda_{1} \left(F_{load} - 2.0 cfr_{kela} v_{L}\right) \cos{\left (\theta \right )} - L \lambda_{3} \left(M_{car} M_{load} g \sin{\left (\theta \right )} + 2.0 M_{car} cfr_{load} v_{vnu} \cos{\left (\theta \right )} + M_{load} \left(F_{car} + 2 F_{load} \sin{\left (\theta \right )} - 2.0 cfr_{car} v_{vnu} - 4.0 cfr_{kela} v_{L} \sin{\left (\theta \right )}\right) \cos{\left (\theta \right )}\right) + \lambda_{2} \left(F_{car} M_{load} \sin{\left (\theta \right )} - 1.0 F_{load} M_{load} \cos{\left (2.0 \theta \right )} - M_{car} M_{load} g \cos{\left (\theta \right )} + 2.0 M_{car} cfr_{load} v_{vnu} \sin{\left (\theta \right )} - 2.0 M_{load} cfr_{car} v_{vnu} \sin{\left (\theta \right )} + 2.0 M_{load} cfr_{kela} v_{L} \cos{\left (2.0 \theta \right )}\right)\right)\\- \frac{1.0}{L^{2} M_{car} M_{load}} \left(1.0 F_{car} M_{load} \lambda_{2} \cos{\left (\theta \right )} + 0.5 F_{load} M_{load} \lambda_{2} \sin{\left (2.0 \theta \right )} + 1.0 L^{2} M_{car} M_{load} \lambda_{3} \omega^{2} + 1.0 M_{car} M_{load} g \lambda_{2} \sin{\left (\theta \right )} + 2.0 M_{car} M_{load} \lambda_{2} \omega v_{L} + 2.0 M_{car} cfr_{load} \lambda_{2} v_{vnu} \cos{\left (\theta \right )} - 2.0 M_{load} cfr_{car} \lambda_{2} v_{vnu} \cos{\left (\theta \right )} - 1.0 M_{load} cfr_{kela} \lambda_{2} v_{L} \sin{\left (2.0 \theta \right )}\right)\end{matrix}\right]$$

Optimiohjauksen numeerinen ratkaiseminen vaatii alkuarvauksen ratkaisulle ajan funktiona. Alkuarvauksen ei helpoissa tapauksissa tarvitse olla kovin hyvä, mutta ratkaisu ei välttämättä lähde konvergoimaan, jos alkuarvaus on esimerkiksi pelkkiä nollia.

Laitan alkuarvaukseksi vakionopeuden vaunulle ja vaijerin pituudelle, minkä perusteella vaunun paikan ja vaijerin pituuden alkuarvaus on helppo laskea.

Liittotiloilla laitoin alkuarvauksen summamutikassa. Arvelin, että joku nollasta poikkeava olisi hyvä.

In [63]:
guess = [(xf_car - x0_car) / T_f, 0.0, (L_f - L_0) / T_f,
         x0_car + (xf_car - x0_car) / T_f * t, 0.0,
         L_0 + (L_f - L_0) * t / T_f, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0]
guess
Out[63]:
$$\left [ \frac{1}{T_{f}} \left(- x_{0 car} + xf_{car}\right), \quad 0.0, \quad \frac{1}{T_{f}} \left(- L_{0} + L_{f}\right), \quad x_{0 car} + \frac{t}{T_{f}} \left(- x_{0 car} + xf_{car}\right), \quad 0.0, \quad L_{0} + \frac{t}{T_{f}} \left(- L_{0} + L_{f}\right), \quad 1.0, \quad -1.0, \quad 1.0, \quad -1.0, \quad -1.0, \quad 1.0\right ]$$

Yhdistetään järjestelmän tilamuuttujat ja liittotilamuuttujat vektoriksi z0 ja niiden aikaderivaatat vektoriksi dz0.

Tässä vaiheessa niiden lausekkeisiin voi sijoittaa optimaalisen ohjauksen lausekkeet kommenteissa esitetyllä tavalla tai voi generoida optimaaliset ohjaukset palauttava funktio osaksi numeerista ratkaisua.

Jälkimmäinen vaihtoehto antaa mahdollisuuden asettaa ohjauksilla helposti ylä- ja alarajat.

In [64]:
z0 = sp.Matrix(yy).col_join(yy_lamda).subs(
    [(F_car, uopt[F_car]), (F_load, uopt[F_load])])
dz0 = sp.Matrix(dx).col_join(dlt).subs(
    [(F_car, uopt[F_car]), (F_load, uopt[F_load])])
z0 = sp.Matrix(yy).col_join(yy_lamda) #.subs([(F_car,uopt[F_car]),(F_load,uopt[F_load])])
dz0 = sp.Matrix(dx).col_join(dlt) #.subs([(F_car,uopt[F_car]),(F_load,uopt[F_load])])

Mikäli ohjauksen lausekkeet sijoittaa yhtälöihin jo tässä vaiheessa, numeerinen algoritmi voi käyttää hyödyksi esimerkiksi tilayhtälön jacobin matriisia.

(Paitsi, että tulos ei kelvannut bvp_solver algoritmille. Se arvelee koodiani virheelliseksi, koska sen numeerisesti laskema arvo ei aivan täsmää koodini tuottaman kanssa.)

jacdz = [[0 for i in range(0,4*sysdim)] for j in range(0,4*sysdim)]
for i in range(0,4*sysdim):
    for j in range(0,4*sysdim):
        jacdz[i][j] = sp.diff(dz[i],z[j])
jacdz

Luodaan lista, jonka avulla voidaan korvata tilamuuttujat vektorilla zz. Koska sympy-lausekkeisiin on vaikea ujuttaa muotoa zz[i] olevia termejä, käytän funktio-muotoa zz(i), jonka korvaan muodolla zz[i] koodin generoinnin yhteydessä.

In [65]:
zz_subs = [(Var, "zz(" + str(i) + ")") for i, Var in enumerate(z0)]
par_subs.extend(zz_subs)
par_subs
Out[65]:
[('sysdim', 'pp(sysdim)'),
 ('udim', 'pp(udim)'),
 ('g', 'pp(g)'),
 ('M_car', 'pp(M_car)'),
 ('M_load', 'pp(M_load)'),
 ('T_f', 'pp(T_f)'),
 ('x0_car', 'pp(x0_car)'),
 ('xf_car', 'pp(xf_car)'),
 ('L_0', 'pp(L_0)'),
 ('L_f', 'pp(L_f)'),
 ('cfr_car', 'pp(cfr_car)'),
 ('cfr_load', 'pp(cfr_load)'),
 ('cfr_kela', 'pp(cfr_kela)'),
 ('C_W', 'pp(C_W)'),
 ('C_F1', 'pp(C_F1)'),
 ('C_F2', 'pp(C_F2)'),
 ('u1lim', 'pp(u1lim)'),
 ('u2lim', 'pp(u2lim)'),
 (v_vnu, 'zz(0)'),
 (omega, 'zz(1)'),
 (v_L, 'zz(2)'),
 (x_car, 'zz(3)'),
 (theta, 'zz(4)'),
 (L, 'zz(5)'),
 (lamda_1, 'zz(6)'),
 (lamda_2, 'zz(7)'),
 (lamda_3, 'zz(8)'),
 (lamda_4, 'zz(9)'),
 (lamda_5, 'zz(10)'),
 (lamda_6, 'zz(11)')]

Kun en ehtinyt muuta ratkaisua hakea, muutan sympy.Matrix-muotoiset matriisit python-listoiksi alkio alkiolta sijoittaen.

Hyvällä onnella simplify sieventää lausekkeita niin, että koodia tulee vähemmän.

In [66]:
z = [sp.simplify(z0i) for z0i in z0]
dz = [sp.simplify(dz0i) for dz0i in dz0]

Koodin generointi

Ratkaisen kahden pisteen reuna-arvotehtävän scikits-paketista löytyvällä bvp_sp.solver-algoritmilla, jota alla generoitava ohjelma kutsuu. scikits dokumentaatiota

Installointi: pip install scikits.bvp_sp.solver

Joidenkin lausekkeiden ympärille olen laittanut "ylimääräiset" sulkeet, jotta python paremmin ymmärtäisi jatkorivejä, joita ne tuottavat.

Generoidaan moduli, jossa on järjestelmän parametrit sisältävä luokkamäärittely Cl_par ja luodaan sitä vastaava olio par. Muissa moduleissa on komento from opt_params import par

Toimii, vaikkei ehkä pitäisi. Olion par importointi ei taida olla pythonin sääntöjen mukaista. Selvittelen.

Koodin generointi write-käskyillä on kömpelöä, mutta ei niin kömpelöä, että olisin jaksanut ryhtyä kehittämään jotain elegantimpaa.

In [67]:
with open('opt_params.py', 'w') as fil:
    fil.write("class Cl_par:\n")
    fil.write("    def __init__(self):\n")
    for (Id, Val) in params:
        s = "        self." + Id + " = " + str(Val) + "\n"
        fil.write(s)
    fil.write("\n")

    fil.write("\npar = Cl_par()\n")
In [68]:
with open('opt_dfuns.py', 'w') as fil:
    fil.write("# -*- coding: utf-8 -*-\n")
    fil.write("#import scikits.bvp_sp.solver\n")
    fil.write("#from scipy.optimize import minimize\n")
    fil.write("import numpy\n")
    fil.write("from math import sin, cos, sqrt \n")
    fil.write("from opt_params import par \n")

    fil.write(" \n")

    fil.write("def uopt(zz):\n")
    fil.write("    F_car = (")
    wr_expr(fil, uopt[F_car], "", par_subs)
    fil.write(")\n    F_load = (")
    wr_expr(fil, uopt[F_load], "", par_subs)
    fil.write(")\n")
    # fil.write(")\n    F_0 = (")
    # wr_expr(fil,F_0,"", par_subs)
    # s = ")\n    F_car = min(par.u1lim,max(-par.u1lim,F_car))\n"
    # fil.write(s)
    # s = "\n    F_load = min(F_0+par.u2lim,max(F_0-par.u2lim,F_load))\n"
    # fil.write(s)
    fil.write("    return [F_car, F_load]\n")
    fil.write(" \n")

    v_krm = v_krm.subs(dx_subs)
    fil.write("def v_kuorma(zz):\n")
    fil.write("    v_kuorma = (")
    wr_expr(fil, v_krm, "", par_subs)
    fil.write(")\n")
    fil.write("    return v_kuorma\n")
    fil.write(" \n")

    fil.write("# Muuttujien derivaatat\n")
    fil.write("def odefun(t,zz):\n")
    fil.write("    [F_car, F_load] = uopt(zz)\n")
    fil.write("    dzdt = ")
    wr_expr(fil, dz, " ", par_subs)
    fil.write("\n    return numpy.array(dzdt)\n")
    fil.write(" \n")

    # Allaoleva antaa bvp_solverin mielestä väärän tuloksen, joten sitä ei kannata generoida
    # fil.write("def gradodefun(t,zz):\n")
    # fil.write("    [F_car, F_load] = uopt(zz)\n")
    # fil.write("    dxdf = ")
    # wr_expr(fil,jacdz,"    ", par_subs)
    # fil.write("\n    return numpy.array(dxdf)\n")
    # fil.write(" \n")

    fil.write("# Reuna-arvot. Palauttaa poikkeaman halutusta.\n")
    fil.write("def bcfun(zz0,zz):\n")
    fil.write("    yy_0 = ")
    wr_expr(fil, yy_0, "", par_subs)
    fil.write("\n    yy_f = ")
    wr_expr(fil, yy_f, "", par_subs)
    fil.write("\n")
    #    fil.write("\n    lamda_f = ")
    #    wr_expr(fil,lamda_f,"", par_subs)
    fil.write("\n    bc0 = [zz0[i]-yy_0[i] for i in range(0,2*par.sysdim)]\n")
    fil.write("\n    bcf = [zz[i]-yy_f[i] for i in range(0,2*par.sysdim)]\n")
    #    fil.write("    bcf = [zz[i+2*par.sysdim]-lamda_f[i] for i in range(0,2*par.sysdim)]\n")
    fil.write("    return(numpy.array(bc0),numpy.array(bcf))\n")
    fil.write(" \n")

    fil.write("# Alkuarvaus ajan funktiona\n")
    fil.write("def init_guess(t):\n")
    fil.write("    guess1 = "),
    wr_expr(fil, guess, "", par_subs)
    fil.write("\n    return guess1\n")
    fil.write(" \n")

    xy_load_lst = [xy.subs(dx_subs) for xy in xy_load]

    fil.write("def xy_load(zz):\n")
    fil.write("    xy_load = ")
    wr_expr(fil, xy_load_lst, " ", par_subs)
    fil.write("\n    return xy_load \n")
    fil.write(" \n")

Esimerkki siitä, miten kopioida tiedoston sisältö toiseen tiedostoon.

with open('opt_piirto.py', 'w') as fil:
    with open('opt_piirto_nogen.py', 'r') as fin:
        fil.write(fin.read())

Hajahuomioita

Käsivoimin yhtälöitä pyörittäville ja koodia kirjoittaville edelläesitetty voi näyttää hämmentävältä ja epäilyttävältäkin. Ohjelma kuitenkin toimii hienosti.

Ylläesitetyn tapaista yhtälöiden pyörittelyä ja koodin generointia voi tehdä myös Maximalla

Matemaattisesta tehtävästä voi selvittää, onko sillä analyyttinen ratkaisu vai ei. Numeerisen ratkaisun löytymisestä ei voi olla varma. Esimerkiksi kahden pisteen reuna-arvotehtävien ratkaisemissa ratkaisun löytyminen saattaa riippua alkuarvauksesta ja tässä tehtävässä optimointikriteerissä käytetyistä kertoimista ja itse järjestelmän parametreistä.

Jos ratkaisu löytyykin, on vaikea tietää, onko se globaali vai lokaali optimi. Jos ratkaisu näyttää hyvältä ja toimii, se saattaa riittää kertaluonteisissa tehtävissä. Jos algoritmin tehtävä olisi esimerkiksi ohjata laiva automaattisesti laituriin, sen pitää antaa joka kerta vähintäänkin kelvollinen ratkaisu.

Tuloksia en isommin tarkastele, koska tarkoitukseni ei ollut tutkia siltanosturin dynamiikkaa vaan kokeilla ja esitellä dynamiikan tutkimisessa hyödyllisiä menetelmiä.

Demovideolla näemme, että optimiohjaus siirtää kuorman siististä paikasta toiseen. Se on melko vahva osoitus siitä, että optimointialgoritmissa ei ole virheitä. Sen perusteella, että simulaattori toimii luonnollisen tuntuisesti, voi arvella, että tilayhtälöiden johtaminen onnistui.

Käytännössä etukäteen laskettu optimiohjaus ei toimisi noin hyvin, koska täydellisen tarkkaa mallia ei voi laatia kuin aivan yksinkertaisille järjestelmille. Lisäksi järjestelmään voi vaikuttaa ulkopuolisia häiriöitä, kuten tuulenpuuskat, joita ei voi ennakoida.

Optimiohjauksen laskemista käytetään kuitenkin hyväksi malliprediktiivisessä säädössä, jossa optimiohjaus lasketaan määrävälein uudestaan lähtien järjestelmän sen hetkisestä tilasta. Järjestelmän kulloinenkin tila määritetään mittaustietoon ja järjestelmän malliin perustuvalla tilaestimaattorilla.

Kotisivuni

Kotisivulleni