Kurs:Mathematische Modellbildung/Themen/2025-26 Wintersemester/Thema 3/Octave-Skript
Octave-Skript mit implementierten Kommentaren
[Bearbeiten]## Anzahl der betrachteten Produkte/Variablen
m=5
## Anzahl der betrachteten Nährstoffe/(Un-)Gleichungen
n=3
## Hilfsvektor der Dimension (m+n)*1, erste (m-1) Einträge =1, Rest =0
V=[1;1;1;1;0;0;0;0]
##(Un-)Gleichungen
# Anzahl der Produkte >=0
M=eye(m)
b=zeros(m,1)
# Nährstoffe der Produkte
f_1=[8.2,14.5,18,5.2,12]
f_2=[78,0.5,15,6.4,71]
f_3=[0.5,8,9.6,1.5,1.4]
# Nährstoffgrenzwerte
ew_1=57
kh_2=287.5
fe_3=70
# M*x>=b
M=[M;f_1;f_2;f_3]
b=[b;ew_1;kh_2;fe_3]
## Kosten der Produkte
c_1=0.299
c_2=0.5725
c_3=1.063
c_4=0.158
c_5=0.138
## Kostenfunktion
c=@(x_1,x_2,x_3,x_4,x_5) c_1*x_1+c_2*x_2+c_3*x_3+c_4*x_4+c_5*x_5
## Hilfsvariable, um zu sehen, ob bessere Lösung gefunden
v1=0
## Hilfsvariable, um zu sehen, ob alle anderen Ecken betrachtet wurden
v2=0
## Hilfsvariable, um Ecke mit anderen Ecken zu vergleichen
v3=0
## Hilfsvariable, um Matrix auf eine andere Ecke zu fixieren
v4=1
## Hilfsvariable, die die zu rotierenden Zeilen der aktuell besten Lösung angibt (1,2,...,(m-1),v5)
v5=0
## mögliche weitere Pfade
Pfadf=[]
## alle Ecken mit gleichguten Kosten zu aktueller bester Lösung
E=[]
for i=m:(m+n) ## erste m-1 Gleichungen festhalten, m-te Gleichung wechseln
H=[M(1,:);M(2,:);M(3,:);M(4,:);M(i,:)] ## Hilfsmatrix
h=[b(1);b(2);b(3);b(4);b(i)] ## Hilfsergebnis
x=H\h ## GLS lösen
if all(M*x>=b-0.001) ## erfüllt x die Ungleichungen? (Rundungsfehler)
Lsg=x ## erste Ecke
cLsg=c(x(1),x(2),x(3),x(4),x(5)) ## Kosten der ersten Ecke
M([1,2,3,4,i],:)=M([i,1,2,3,4],:) ## verschiebt Zeile i nach oben
b([1,2,3,4,i])=b([i,1,2,3,4]) ## verschiebt ebenso
rotationen=1 ## erste Rotation gerade durchgeführt
E=x' ## Ecke zu gefundenen Ecken hinzufügen
v5=i
break ## beendet die for-Schleife
endif
endfor
while v2==0
while rotationen<m ## solange noch nicht ganz durchrotiert
for i=m:(m+n) ## erste m-1 Gleichungen fest
v1=0
H=[M(1,:);M(2,:);M(3,:);M(4,:);M(i,:)] ## Hilfsmatrix
h=[b(1);b(2);b(3);b(4);b(i)] ## Hilfsergebnis
x=H\h ## GLS lösen
if all(M*x>=b-0.001) ## erfüllt x die Ungleichungen? (Rundungsfehler)
if c(x(1),x(2),x(3),x(4),x(5))<cLsg ## ist x eine besssere Lösung?
Lsg=x ## bessere Lösung abspeichern
cLsg=c(x(1),x(2),x(3),x(4),x(5)) ## Kosten der besseren Lösung
M([1,2,3,4,i],:)=M([i,1,2,3,4],:) ## Matrix durchrotieren
b([1,2,3,4,i])=b([i,1,2,3,4]) ## entsprechend b rotieren
rotationen=1 ## Rotation zurücksetzen (werden nach for-Schleife um 1 erhöht)
v1=1
Pfadf=[] ## weitere Pfade können ignoriert werden
v5=i
E=x' ## Ecke zu gefundenen Ecken hinzufügen
break ## beendet die for-Schleife
elseif c(x(1),x(2),x(3),x(4),x(5))<cLsg+0.001 ## ist x eine genausogute Lösung (mit Rundungsfehler)?
for k=1:rows(E) ## überprüfe, ob Ecke bereits gefunden
if all(abs(x' - E(k,:))<0.001) ## Ecke wurde gefunden, mit Rundeungsfehler
v3=1
endif
endfor
if v3==0 ## Ecke wurde noch nicht gefunden
E=[E;x'] ## Ecke zu gefundenen Ecken hinzufügen
Pfadf=[Pfadf,[V+eye(m+n)(:,i)]] ## möglichen Pfad für später in Spaltenvektor speichern
## erzeugt Standard-Einheitsvektor
endif
v3=0
endif
endif
endfor
if v1==0 ## keine bessere Ecke wurde gefunden --> Rotieren
M([1,2,3,4,v5],:)=M([v5,1,2,3,4],:) ## Matrix rotieren
b([1,2,3,4,v5])=b([v5,1,2,3,4]) ## Vektor rotieren
rotationen=rotationen+1
if !isempty(Pfadf) ## damit kein Fehler in Octave
Pfadf([1,2,3,4,v5],:)=Pfadf([v5,1,2,3,4],:) ## Pfadgleichungen rotieren
endif
endif
endwhile
if isempty(Pfadf) ## keine weiteren Ecken mit den selben Kosten
v2=1
else
rotationen=0
for i=1:(m+n) ## Matrix anpassen, damit Algorithmus für neue Ecke starten kann
if !isempty(Pfadf) ## damit kein Fehler in Octave
if Pfadf(i,1)==1 ## nur verwendete Gleichungen werden betrachtet
M([v4,i],:)=M([i,v4],:) ## i-te Zeile kommt in v4-te Zeile
Pfadf([v4,i],:)=Pfadf([i,v4],:) ## entsprechend Pfadgleichungen anpassen
v4=v4+1
endif
endif
endfor
v4=1
if !isempty(Pfadf) ## damit kein Fehler in Octave
Pfadf(:,1)=[] ## Ecke/Pfadgleichungen entfernen, da diese jetzt betrachtet werden
endif
v5=m ## Festlegen der zu rotierenden Zeilen 1,2,3,...,m
endif
endwhile
## Ausgabe der besten Lösung
Lsg
cLsg
Zum Debuggen des Programmes wurde Chat AI|Academic Cloud [1] verwendet.
Anpassungen des Octave-Skripts
[Bearbeiten]Das Octave-Skript kann in wenigen Schritten per Hand angepasst werden, falls eine andere Anzahl an Produkten (Variablen) oder Nährstoffen (Gleichungen) betrachtet wird:
| Zeile | zu verändernder Text im Skript |
|---|---|
| 2 | Anzahl der Produkte verändern |
| 4 | Anzahl der Nährstoffe verändern |
| 7 | Vektor mit (m-1) 1en und (n+1) 0en |
| 14-16 | Nährstoffe der Produkte anpassen (erster Eintrag Produkt x_1, …), Zeilenvektoren hinzufügen |
| 18-20 | Nährstoffgrenzwerte anpassen, Grenzwerte hinzufügen |
| 22-23 | Gleichungen bzw. Nährstoffgrenzwerte in Matrix hinzufügen |
| 26-30 | Kosten der Produkte anpassen, Kosten hinzufügen |
| 32 | c=@(x_1, x_2, …, x_m) c_1*x_1+c_2*x_2+…+c_m*x_m |
| 51,70 | Matrix H=[M(1,:), M(2,:), …, M((m-1),:), M(i,:)] |
| 52,71 | Vektor h=[b(1), b(2), …, b(m-1), b(i)] |
| 56,74,76,85 | c(x(1), x(2), …, x(m)) |
| 57,77 | M([1 ,2 ,… ,(m-1) ,i],:)=M([i ,1 ,2 ,… ,(m-1)],:) |
| 58,78 | b([1 ,2 ,… ,(m-1) ,i])= b([i, 1 ,2 ,… ,(m-1)]) |
| 101 | M([1 , 2, …, (m-1), v5],:)=M([v5, 1 , 2, …, (m-1)],:) |
| 102 | b([1 , 2, …, (m-1), v5])=b([v5, 1 , 2, …, (m-1)]) |
| 105 | Pfadf([1 , 2, …, (m-1), v5],:)=Pfadf([v5, 1 , 2, …, (m-1)],:) |
Zum Anzeigen der Zeilennummern im Programm, muss dieses lediglich in ein Octaave-Skript kopiert werden.
Zurück zur Modellierung: Ein optimaler Ernährungsplan