Inleiding Programmeren + R

for-loop voorbeeld

» Start

Een voedselregulatie model


Veel technologische feedbacksystemen worden gecontroleerd door een mechanisme dat een proces start of stopt. De thermostaat van een verwarming is hier een goed voorbeeld van. Wanneer het in een kamer te koud wordt dan zorgt de thermostaat ervoor dat de verwarming aanslaat. Wanneer het in een kamer te warm wordt dan zorgt de thermostaat dat de verwarming uit gaat.


Een biologisch voorbeeld van een dergelijk controlesysteem is het voedingsgedrag van ratten. In een sterk vereenvoudigd model is in het voedsel systeem een regelmechanisme ingebouwd dat een 'switch' op aan zet als de rat gaat eten en de switch op uit zet als de rat niet eet.

Van dit probleem is onderstaand diagram gemaakt. 


feedback.png

 

De darm dient in dit model als voedsel opslagplaats. Via de darmwand vindt energie toevoer naar het lichaam plaats. De snelheid van energie transport naar het lichaam is evenredig met het oppervlak van de darmwand. In dit model nemen we aan dat dit evenredig is met de wortel van de darm inhoud:


M = km √D


Hier is M de snelheid van energie overdracht van de darm naar het lichaam, km is de opnamecoëfficiënt en D is de hoeveelheid energie van het voedsel in de darmen. km kan gedurende een etmaal variëren. Omdat de rat een nachtdier is zal de coëfficiënt 's nachts hoger zijn dan overdag. Opname van voedsel vindt snel plaats, terwijl energieverbruik langzaam gaat. De snelheid van voedselopname, F, wordt berekend door:


F = eI


waarbij I de eetsnelheid is en e de feedback coëfficiënt (het regel mechanisme). Als de rat eet, is e gelijk aan 1. De darm stroomt dan met snelheid I vol. Als de rat niet eet is e gelijk aan 0.

In ons model wordt het regelmechanisme aan- en uitgezet door de hoeveelheid energie M. Als M onder een bepaald niveau komt, Ml, dan zal de rat gaan eten. Als M boven een maximum komt, Mh, dan stopt de rat met eten en wordt e op 0 gezet.

In dit model wordt aangenomen dat het gewicht van de rat constant blijft. De energie kan niet worden opgeslagen, alleen energie uit de darmen wordt gebruikt. Voedsel is ten alle tijden aanwezig. 

De uiteindelijke formule voor het systeem is:


D(t) = D(t-1) + F - M


We maken een algoritme waarin we 24 uur (1440 minuten) van de rat simuleren. We gebruiken tijdstapjes van 5 minuten.


Gegeven zijn de volgende startwaarden:

Om de 8 uur verandert de km van de rat (variërend van 1 naar 2 naar 3), afhankelijk van ochtend, dag, of nacht

Variabelen en hun startwaarden: 

Ml = 18 cal/min

Mh = 60 cal/min

I = 500 cal/min

D = 4000 cal

Begintijd is 7 uur 's ochtends


Dit levert ons het volgende algoritme:

    

    voedselregulatie <- function() {

        I <- 500

        Ml <- 18

        Mh <- 60

        D <- 4000

        km <- 1                     #km op 'ochtendstand' 

        e <- 0

        M <- km * sqrt(D)

        for (tijd in 1:288) {       #Doe 288 tijdstappen 

            if (M < Ml) {

                e <- 1              # zet eetsysteem aan

                }

            if (M > Mh) {

                e <- 0              # zet eetsysteem uit

                }

            F <- e * I              # snelheid van darm-vullen 

            D <- D + F - M          # bereken darminhoud 

            M <- km * sqrt(D)       # energie overdracht 

            cat(c(F, D, M),"\n")    # schrijf diverse model parameters

            if (tijd == 96) {       #om de 8 uur veranderd km 

                km <- 2             # km op 'dagstand' 

                }

            if (tijd == 192) { 

                km <- 3             # km op 'nachtstand'

                }

        }                           # einde for lus

    }                               # einde functie definitie

    

‘Prachtig hoor’, zul je zeggen, ‘maar hoe kom ik hier nu aan als ik het algoritme zelf moet schrijven ?  Geef eens wat richtlijnen ! ‘


Hoe bouwen we het algoritme op ?


De beginregel is het eenvoudigst. Dat is de programma-definitie regel. Op die regel wordt de naam van het programma vastgelegd en, als daarom is gevraagd, ook de mogelijkheid om input (= parameters) mee te geven. In dit geval wordt er niet gevraagd om parameters over te kunnen dragen. Dat betekent dat onze programma definitie regel er als volgt uit moet zien: 

voedselregulatie <- function() {

De volgende stap is de initialisatie-fase, en bestaat er uit dat we ons afvragen hoeveel variabelen er in het probleem dat we moeten oplossen een rol spelen, van welk type ('getal' of 'teken') ze zijn, als ze het type 'getal' hebben welke vorm (dimensie) ze moeten hebben (scalair, vector, of matrix: zie gegevensstructuren), en welke startwaarde de variabelen bij het begin van het programma moeten krijgen.


In dit geval wordt het aantal variabelen bepaald door de formule waarmee we de toestand van het metabolisme per tijdstap uitrekenen. In die formule spelen een zevental variabelen een rol. 

In de eerste plaats de eetsnelheid I, vervolgens de snelheid van energie overdracht M, het minimum- en maximum niveau, Ml en Mh, daarvan, de opname coëfficient km, de hoeveelheid energie D van het voedsel in de darm, en tot slot de feedback coëfficient e. Daarmee is het lijstje van variabelen kompleet. Alle zeven hebben ze aan het begin een startwaarde nodig; ze moeten worden geinitialiseerd anders kunnen we er verder op in het programma niet mee rekenen. Verder moeten we er rekening mee houden dat de waarde van de variabele M niet direct bekend is maar afhankelijk is van de waarde van km en D. M kan dus pas worden uitgerekend nadat km en D hun waarde hebben gekregen ! Het initialisatie blok komt er dan als volgt uit te zien:

 

I  <-  500

Ml  <-  18

Mh  <-  60

D  <-  4000

km  <-  1 #km op 'ochtendstand' 

e  <-  0

M  <-  km * sqrt(D)


OK, de voorbereidende stappen zijn achter de rug en we zijn nu aan het echte werk toe. Wat is het probleem ? We moeten voor een vaststaande periode (24 uur met 1440 minuten) in een bepaalde regelmaat (om de 5 minuten) volgens bepaalde formules de waarden van variabelen bepalen. De kern van het programma moet dus een lus bevatten, want we moeten om de 5 minuten steeds dezelfde formules toepassen. En aangezien we precies weten hoeveel blokjes van 5 minuten er in 24 uur gaan (288) ligt het soort lus ook vast, nl de voor-lus ! De structuur van een voorlus is als volgt:


for (tijdstip in 1:totaalAantalTijdstippen) {

opdracht 1

opdracht 2

opdracht 3

...

opdracht n

}


In R wordt dat:

for (tijd in 1: 288 ) { #Doe 288 tijdstappen 

We weten ook wat we in moeten vullen voor het blokje opdrachten, namelijk de tot geldige R expressies omgeschreven formules ! Dat ziet er als volgt uit:

F  <-  e * I # snelheid van darm-vullen 

D  <-  D + F - M # bereken darminhoud 

M  <-  km * sqrt(D) # energie overdracht 

Eenmaal uitgerekend moet de gebruiker van het programma worden verteld wat nu de waarden van de verschillende variabelen is geworden. Dat kunnen we doen door die waarden naar het scherm te schrijven met een SCHRIJF opdracht: 

cat(c(F, D, M),"\n") # schrijf diverse model parameters

Daarmee is de kous nog niet af ! De waarde van een aantal van de variabelen is namelijk niet afhankelijk van de formule maar van het tijdstip van de dag (km) of van een bepaalde limiet overschrijding (e). Dat betekent dat we nog een aantal voorwaardelijke opdrachten moeten formuleren waarin die wisselende omstandigheden tot uitdrukking worden gebracht. De eerste serie voorwaarden betreffen de snelheid van energie overdracht. Zodra we voor de eerste keer M hebben uitgerekend (in het initialisatie blok) moeten we controleren of de minimum waarde Ml dan wel de maximum waarde Mh van M is overschreden. In die gevallen moeten we namelijk de waarde van de feedback coëfficient e wijzigen. Direct na het begin van de voorlus, maar nog voor het formule blok, komen deze voorwaarden te staan:

if  (M < Ml)  {

e  <-  1 # zet eetsysteem aan

}

if (M > Mh)  {

e <-  0 # zet eetsysteem uit

}

De andere variabele die voorwaardelijk verandert is km, de opnamecoëfficient. Deze verandert afhankelijk van het tijdstip. We weten dat we de simulatie om 7 uur ‘s ochtends beginnen, en dus weten we de begin toestand van km, nl de ochtendstand. Maar km moet veranderen zodra we 8 uur verder zijn, en weer een keer na nog eens 8 uur. Die voorwaardelijke veranderingen in de waarde van km komen na het formule blok. (maakt dat wat uit, of we km voor of na het formule blok controleren ?):

if (tijd == 96) { #om de 8 uur veranderd km 

km  <-  2 # km op 'dagstand' 

    }

if (tijd == 192) { 

km  <-  3 # km op 'nachtstand'

}

Hiermee zijn de bouwstenen van het algoritme kompleet. Nu alles nog in de goede volgorde zetten.