Hallo zusammen,
„i produly present“:
- Machine Learning Script zur Vorhersage der PV Daten mit einem Neuronalen Netz.
Da ich noch keine Langzeit-Daten verfügbar habe, könnt Ihr das gerne mal testen und Feedback geben nach einigen Tagen.
Ich hab alles in ein Scriptfile gepackt, so dass nur das Skript angelegt werden muss. Der Rest sollte nach Konfiguration der Parameter automatsich funktionieren.
Bitte im Kopf des Skriptes die Parameter setzen.
<?php
/* -------------------------------------------------------------
Machine Learning PV Forecast
(c) 2022/09 stele99
Skript baut ein Neuronales Netzwerk auf um PV Ertragsdaten
vorherzusagen. Dabei ist nur die Groesse, der Standort
anzugeben und die naechste DWD Messtation anzugeben.
Das NN wird mit dem Sonnenstand und der Globalstrahlung gefuettert
und mit dem echten Ertrag als Zielwert trainiert.
Damit das Training durchgefuehrt werden kann, werden 2
Archivariablen angelegt, die stuendlich aktualisiert werden.
Auf Basis dieser Daten kann dann der Forecast der PV Anlage
(unabh. von Verschattung / Ausrichtung) ermittelt werden.
Der Forecast wird in der Variable "JSON PV Forecast" abgelegt.
Dies ist ein JSON Array mit den Forecast und WEtterdaten
Der stuendliche Forecast liegt in dem Feld pv_estimate
Das Script sollte erstmal 1-3 Tage Daten sammeln.
Last Udpate: 2022-09-14 19:54 Uhr
---------------------------------------------------------------*/
/* ### Einstellungen ########################################### */
$PV_MAX = 8000; // Watt Peak der PV Anlage
$PV_CURRENT = 39680; // Variable mit aktueller Leistung der PV Anlage - sollte moeglichst realtime sein.
$PV_LAT = 48.024873309; // Breitengrad der Anlage
$PV_LON = 8.929772606; // Laengengrad der Anlage
$FC_STATION = 10818; // Stationsid der DWD Vorhersage https://wettwarn.de/mosmix/mosmix.html
$ID_ARCHIVE = 26472;// ID des Archive Handlers
$GSMAX = 1800; // Maximale Globalstrahlung
$NN_NAME = "pv_garage"; // Name für neuronales Netz (=Dateiname für Persistenz), kann frei geaendert werden.keine sonderzeichen
$FC_VARIABLE = true; // Erzeuge Variable fuer forecast in Watt fuer heute und morgen
$FC_COMPARE = true; //true wenn PV-Forecast in Variable gespeichert werden sollen fuer spaeteren Vergleich (SOLL /IST)
$FC_CHART = true; // Generiere Vorhersagechart fr den aktuellen Tag und Morgen
$NN_TRAIN = 1; // 1 = trainieren aktiv, 0 kein Training mehr, -1 Training löschen
$NN_TRAIN_MAX = 365; // Anzahl der Tage für Training // Achtun CPU Bedarf steigt!
/* #### weitere CONFIG Daten (intern) ################################### */
$FC_CHACHEFILE = dirname(__FILE__) . "/cache/".$_IPS['SELF']."-".$FC_STATION.".cache";
$FC_CHACHEAGE = 3600 * 3 ; // WetterCache 3 Stunden // Frage DWD nur alle 3 Stunden ab.
$FC_XMLFILE = dirname(__FILE__) . "/cache/".$_IPS['SELF']."-".$FC_STATION.".xml"; // XML Forecast entzippt
$FC_APIURL = "http://opendata.dwd.de/weather/local_forecasts/mos/MOSMIX_L/single_stations/" . $FC_STATION . "/kml/MOSMIX_L_LATEST_" . $FC_STATION . ".kmz";
/* ### Skript regelmäßig laufen lassen, Ereignis anlegen, falls nicht vorhanden### */
$subIDs = IPS_GetChildrenIDs($_IPS['SELF']);
$evt = false;
foreach($subIDs as $id){
$obj = IPS_GetObject($id);
if($obj["ObjectType"] == 2)$evt = true;
}
if(!$evt){
$eid = IPS_CreateEvent(1);
IPS_SetName($eid,"Stündliche Aktualisierung und Protokollierung");
IPS_SetParent($eid, $_IPS['SELF']);
IPS_SetEventCyclicTimeFrom($eid, 0, 10, 0);
IPS_SetEventCyclic($eid,0,0,0,0,3,1);
IPS_SetEventActive($eid, true);
}
/* ### DWD FORECAST ########################################### */
//cache verzeichnis anlegen, wenn noch nicht da:
if(! is_dir(dirname(__FILE__) . "/cache")){
mkdir(dirname(__FILE__) . "/cache");
}
// Cache der Abfrage
$lastrun = intval(@filemtime($FC_CHACHEFILE));
if (time() - $lastrun > $FC_CHACHEAGE ) {
$response = file_get_contents($FC_APIURL);
file_put_contents($FC_CHACHEFILE, $response);
$zip = new ZipArchive;
$res = $zip->open($FC_CHACHEFILE);
if ($res === TRUE) {
$zc = $zip->statIndex(0);
$zf = $zc["name"];
$zip->extractTo(dirname(__FILE__) . "/cache/", $zf);
$zip->close();
copy(dirname(__FILE__) . "/cache/".$zf, $FC_XMLFILE);
unlink(dirname(__FILE__) . "/cache/".$zf);
} else {
echo 'Fehler, Code:' . $res;
}
}
$xmlstr = file_get_contents($FC_XMLFILE);
// Namespace aufräumen
$xmlstr = str_replace("<kml:", "<", $xmlstr);
$xmlstr = str_replace("</kml:", "</", $xmlstr);
$xmlstr = str_replace("<dwd:", "<", $xmlstr);
$xmlstr = str_replace("</dwd:", "</", $xmlstr);
$xmlstr = str_replace(" dwd:", " ", $xmlstr);
$xml = simplexml_load_string($xmlstr);
$ts = $xml->Document->ExtendedData->ProductDefinition->ForecastTimeSteps->TimeStep;
$fca["info"]["model"]=$xml->Document->ExtendedData->ProductDefinition->ProductID->__toString();
$fca["info"]["generation_time"]=$xml->Document->ExtendedData->ProductDefinition->IssueTime->__toString();
foreach ($ts as $t) {
$fc = ["ts" => strtotime($t),
"tiso" => $t->__toString(),
"day" => date("z", strtotime($t)) - date("z"),
"hour" => date("G", strtotime($t))];
$fca["hourly"][] = $fc;
}
// Daten aus XML lesen (siehe oben verlinktes excel), 1 Parameter Wertekennung von DWD, 2. Parameter name im Json.
getDataKML("Rad1h", "radiation");
getDataKML("RRad1", "radiation_intensity");
getDataKML("Neff", "clouds");
getDataKML("RRL1c", "prec");
getDataKML("RRS1c", "snow");
getDataKML("SunD1", "sun");
getDataKML("T5cm", "temp");
// Tageswerte berechnen
$dayOld = 0;
$t_min = 999;
$t_max = -999;
$t_avg = 0;
$cnt = 0;
$cloud_avg = 0;
$prec = 0;
$snow = 0;
$sun = 0;
foreach($fca["hourly"] as $fc){
$day = date("z", $fc["ts"]) - date("z");
if($day != $dayOld){
$fca["daily"][$dayOld]["ts"] = mktime(0,0,0, date("m",$fc["ts"]), date("d",$fc["ts"]), date("Y",$fc["ts"]) );
$fca["daily"][$dayOld]["txtx"] = date("D d.m.Y", $fca["daily"][$dayOld]["ts"]);
$fca["daily"][$dayOld]["temp_max"] = $t_max;
$t_max = -999;
$fca["daily"][$dayOld]["temp_min"] = $t_min;
$t_min = 999;
$fca["daily"][$dayOld]["temp_avg"] = round($t_avg / $cnt,1);
$t_avg = 0;
$fca["daily"][$dayOld]["cloud_avg"] = round($cloud_avg / $cnt);
$cloud_avg = 0;
$fca["daily"][$dayOld]["prec"] = $prec;
$prec=0;
$fca["daily"][$dayOld]["snow"] = $snow;
$snow=0;
$fca["daily"][$dayOld]["sun"] = round($sun / 60);
$sun=0;
$cnt = 0;
}
$t_min = ($fc["temp"] < $t_min)? $fc["temp"] : $t_min;
$t_max = ($fc["temp"] > $t_max)? $fc["temp"] : $t_max;
$cnt++;
$t_avg += $fc["temp"];
$cloud_avg += $fc["clouds"];
$prec += $fc["prec"];
$snow += $fc["snow"];
$sun += $fc["sun"];
$dayOld = $day;
$tiso = $fc["tiso"];
}
// Variable erzeugen und Daten als JSON String ausgeben.
$json = json_encode($fca);
$id= CreateVariableByName($_IPS['SELF'],"JSON DWD Forecast", 3);
SetValue($id,$json);
/* #### Historie speichern fuer ML Training #################### */
// Wirkleistung
$id = CreateVariableByName($_IPS['SELF'], "PV Wirkleistung", 2);
$v = IPS_GetVariable($id);
AC_SetLoggingStatus($ID_ARCHIVE, $id, true);
AC_SetAggregationType($ID_ARCHIVE, $id,0);
if(date("G",$v['VariableUpdated']) < date("G")-1 || date("z",$v['VariableUpdated']) < date('z')){
// Wirkleistung der vergangenen Stunde ermitteln (avg) und speichern
$hour = date("G") - 1;
$t1 = mktime($hour,0,0,date('m'),date('d'),date('y'));
$t2 = mktime($hour,59,0,date('m'),date('d'),date('y'));
$w = AC_GetAggregatedValues( $ID_ARCHIVE, $PV_CURRENT, 0, $t1, $t2, 10);
if(!is_array($w))die("Bitte Logging für PV aktivieren ($PV_CURRENT)");
$w = $w[0]["Avg"];
SetValue($id,$w);
}
// Globalstrahlung der vergangen Stunde
$id = CreateVariableByName($_IPS['SELF'], "Globalstrahlung", 2);
$v = IPS_GetVariable($id);
AC_SetLoggingStatus($ID_ARCHIVE, $id, true);
AC_SetAggregationType($ID_ARCHIVE, $id,0);
$radieff = 0;
if(date("G",$v['VariableUpdated']) < date("G")-1 || date("z",$v['VariableUpdated']) < date('z')){
// Globalstrahlung der vergangenen Stunde ermitteln aus Forecast Daten
$hour = date("G"); // DWD gibt letzte Stunde an
$fc = json_decode(GetValue(IPS_GetVariableIDByName("JSON DWD Forecast", $_IPS['SELF'])),true);
foreach($fc["hourly"] as $fh){
if(date("z", $fh["ts"]) == date("z") && date("G", $fh["ts"]) == date("G") ){
$radieff = $fh["radiation"] * ($fh["radiation_intensity"]/100); // effektive Strahlung
break;
}
}
SetValue($id,$radieff);
}
### Neuronales Netz Trainieren, wenn genügend Daten vorhanden sind.
$fn = dirname(__FILE__) . "/".$NN_NAME;
$n = new NeuralNetwork(3, 2, 1); // 3 input / 2 Hidden Layer Neuronen und 1 Output Neuron
if(file_exists($fn)) $n->load($fn);
if($NN_TRAIN == -1)unlink($fn);
if(!file_exists($fn) || ( (filemtime($fn) < time()-3600*24) && $NN_TRAIN == 1 ) ){
// Historie zum Trainieren ermitteln
$idGS = IPS_GetVariableIDByName("Globalstrahlung",$_IPS['SELF']);
$idPV = IPS_GetVariableIDByName("PV Wirkleistung",$_IPS['SELF']);
#$idGS = 39052;
#$idPV = 46621;
$d = $NN_TRAIN_MAX;
$d = $d * 24;
$samples=[];
$targets=[];
for($dt = $d; $dt >0; $dt--)
{
$t = time() - ($dt * 3600);
$hour =date("G", $t);
$t1 = mktime($hour,0,0,date('m',$t),date('d',$t),date('y',$t));
$t2 = mktime($hour,59,0,date('m',$t),date('d',$t),date('y',$t));
// Aus dem Archivehandler die letzte Stunde lesen (Ertrag PV)
$w = AC_GetAggregatedValues( $ID_ARCHIVE, $idPV, 0, $t1, $t2, 10 );
if(count($w)>0){
$w = round($w[0]["Avg"]/10,0)*10;
$gs = round(AC_GetAggregatedValues( $ID_ARCHIVE, $idGS, 0, $t1, $t2, 10 )[0]["Avg"]/10,0)*10;
// Azimuth und Elevation zur 'Stichstunde'
$tav = mktime($hour,30,0,date('m',$t),date('d',$t),date('y',$t));
$sunpos = calcSun($tav, $PV_LAT, $PV_LON);
$azimuth = round($sunpos["azimuth"]/10,0)*10;
$elevation = round($sunpos["elevation"],0);
// Trainieren nur für mögliche Zeiträume
if($hour > 5 && $hour < 22){
$samples[] = [$azimuth, $elevation, $gs];
$targets[] = $w;
}
}
}
if(count($targets)< 5 ){echo "Noch zu wenig Trainingsdaten.";return;};
if(count($targets) < 48) echo "Achtung, noch nicht ausreichend Trainingsdaten...\n";
// Normalisierung der Daten:
//1. Targets
foreach($targets as $k => $t){
$targets[$k] = $t / $PV_MAX;
}
//2. Samples
foreach($samples as $k => $s){
$samples[$k][0] = $samples[$k][0] / 360;
$samples[$k][1] = $samples[$k][1] / 90;
$samples[$k][2] = $samples[$k][2] / $GSMAX;
}
$n->setVerbose(false);
$n->setLearningRate(0.05);
$n->setMomentum(0.3);
$max = 4; // Max Train Runs
$i = 0;
for($x=0;$x < count($samples)-1; $x++){
$n->addTestData($samples[$x], [$targets[$x]]);
}
echo "Start Neural Network Training...\n";
while (!($success = $n->train(10000, 0.06)) && ++$i<$max) {
echo " Train Round $i: No success...\n";
}
if ($success) {
$epochs = $n->getEpoch();
echo " Train Success in $epochs training rounds!\n";
}else{
echo " Train not Successful - using fuzzy data\n";
}
$n->save($fn);
} // NN traning
#### PV Forecast mit NN berechnen ##############################
$fc = json_decode(GetValue(IPS_GetVariableIDByName("JSON DWD Forecast", $_IPS['SELF'])),true);
foreach($fc["hourly"] as $k => $fh){
### ACHTUNG PRÜFEN AKTUELLE STUNDE VS LETZTE STUNDE!!!!!!!!!!!!!!!!!
$ts = mktime(date("H",$fh["ts"]), 30, 0, date("m",$fh["ts"]), date("d",$fh["ts"]), date("Y",$fh["ts"]));
$sunpos = calcSun($ts , $PV_LAT, $PV_LON);
$azimuth = round($sunpos["azimuth"]/10,0)*10 / 360;
$elevation = round($sunpos["elevation"],0) / 90;
$radieff = $fh["radiation"] * ($fh["radiation_intensity"]/100) / $GSMAX; // effektive Strahlung
$pv_estimate = round($n->calculate([$azimuth, $elevation, $radieff])[0] * $PV_MAX/10)*10 ;
if($pv_estimate <0 || $fh["hour"]>21 || $fh["hour"] < 6)$pv_estimate = 0;
$fc["hourly"][$k]["pv_estimate"] = $pv_estimate;
# echo date("Y-m-d, H:",$fh["ts"])." ".$pv_estimate."W\n";
}
$id = CreateVariableByName($_IPS['SELF'], "JSON PV Forecast", 3);
$val = json_encode($fc);
SetValue($id,$val);
// Tageswerte in Variable Werte ausgeben
if($FC_VARIABLE){
$id = CreateVariableByName($_IPS['SELF'], "Erwarteter Ertrag Heute", 2, "~Watt");
$wsum = 0;
foreach($fc["hourly"] as $fh){
if(date("z", $fh["ts"]) == date('z'))$wsum += $fh["pv_estimate"];
}
// nur vor 8 Forecast ausgeben
if(date('G')<8)SetValue($id,$wsum);
$id = CreateVariableByName($_IPS['SELF'], "Erwarteter Ertrag Morgen", 2, "~Watt");
$wsum = 0;
foreach($fc["hourly"] as $fh){
if(date("z", $fh["ts"]) == date('z',time()+3600*24))$wsum += $fh["pv_estimate"];
}
SetValue($id,$wsum) ;
}
// Forecast Comapare einrichten, wenn gewnscht
$fn_dayfc = (__FILE__) . ".fc.cache.json";
if($FC_COMPARE && date('G') < 9 || !file_exists($fn_dayfc)){
file_put_contents($fn_dayfc,$val);
}
$fc_old = json_decode(file_get_contents($fn_dayfc),true); // 1. Forecast des Tages nutzen
// Forecast in Variable speichern falls Compare gewnscht
if($FC_COMPARE){
$id = CreateVariableByName($_IPS['SELF'], "PV Forecast", 2);
$v = IPS_GetVariable($id);
AC_SetLoggingStatus($ID_ARCHIVE, $id, true);
AC_SetAggregationType($ID_ARCHIVE, $id,0);
if(date("G",$v['VariableUpdated']) < date("G")-1 ){
// Globalstrahlung der vergangenen Stunde ermitteln aus Forecast Daten
$hour = date("G"); // werte der letzten Stunden (DWD Format)
foreach($fc_old["hourly"] as $fh){
if(date("z", $fh["ts"]) == date("z") && date("G", $fh["ts"]) == date("G") ){
$pv_estimate = $fh["pv_estimate"];
break;
}
}
SetValue($id,$pv_estimate);
}
}
### Forecast charts #############################################
if($FC_CHART){
$idChart[0] = CreateVariableByName($_IPS['SELF'], "PV Forecast Chart heute", 3, "~HTMLBox");
$idChart[1] = CreateVariableByName($_IPS['SELF'], "PV Forecast Chart morgen", 3, "~HTMLBox");
for($day=0; $day<2; $day ++){
// Read Values from Forecast
$html = '<div style="height: 350px; "><link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/charts.css/dist/charts.min.css">
<style>
#pvfc_day {
--primary-axis-color: rgba(255, 255, 255, 0.3);
}
</style>
';
$html .= '<table class="charts-css line multiple show-data-on-hover reverse-datasets show-labels show-primary-axis" id="pvfc_day">';
$html .= " <caption>PV Forecast today</caption>";
$data_th = "";
$data_fc = "";
$szP = 0;
// Read current data from variable and past
$maxWatt = 0;
if($day == 0){
$idPV = IPS_GetVariableIDByName("PV Wirkleistung",$_IPS['SELF']);
$hour = date("G") - 1;
for($ch = 6; $ch < date("G"); $ch++){
$t1 = mktime($ch,0,0,date('m'),date('d'),date('y'));
$t2 = mktime($ch,59,0,date('m'),date('d'),date('y'));
$w = AC_GetAggregatedValues( $ID_ARCHIVE, $PV_CURRENT, 0, $t1, $t2, 10);
$hwA[$ch+1] = round($w[0]["Avg"]);
if($hwA[$ch+1] > $maxWatt) $maxWatt = $hwA[$ch+1];// Scaling rausfinden
}
}
$szIstP = 0;
// Scaling rausfinden
foreach($fc["hourly"] as $fch){
if($fch["pv_estimate"] > $maxWatt) $maxWatt = $fch["pv_estimate"];
}
$szCloudP=0;
$total = 0;
foreach($fc_old["hourly"] as $fch){
if(date("G", $fch["ts"]) > 6 && date("G", $fch["ts"]) < 22 && date("z", $fch["ts"]) == date('z')+$day){
$hour = date("G", $fch["ts"]);
//Forecast
$data_th .= "<th scope='col'>$hour</th>";
$sz = round($fch["pv_estimate"] / $maxWatt,1);
$data_fc .= '<tr>
<th scope="row"> '.($hour - 1).' </th>
<td style="--start: '.$szP.'; --size: '.$sz.'; --color: #D9A505; font-size:10px; padding-bottom:20px; color: #D9A505">'.$fch["pv_estimate"].'</td>
';
$total += $fch["pv_estimate"];
if(isset($hwA[$hour])&& $day ==0){
$szIst = round($hwA[$hour] / $maxWatt,1);
$data_fc .= '<td style="--start: '.$szIstP.'; --size: '.$szIst.'; --color: #FFFFFF; font-size:10px; padding-bottom:30px;color: #FFFFFF;">'.$hwA[$hour].'</td>';
$szIstP = $szIst;
}
#echo "$hour : ".$hwA[$hour]. " -> " .$fch["pv_estimate"]."\n";
$data_fc .='</tr>';
$szP = $sz;
}
}
$html .= "<thead><tr>$data_th</tr></thead>";
$html .= $data_fc;
$html .= "</table>
<div style='position:relative;top: -350px;'>Gesamtertrag: ".round($total/1000,1)."kw</div>
</div>";
setValue($idChart[$day], $html);
} // Day
}
### HELPER FUNCTIONS ############################################
// Extraktion von Daten us der DWD KML Datei
function getDataKML($idstr, $idtxt)
{
global $xml;
global $fca;
$gs = $xml->Document->Placemark->ExtendedData;
foreach ($gs->Forecast as $g) {
$id = $g->attributes()["elementName"][0]->__toString();
if ($id == $idstr) {
$val = $g->value->__toString();
$valA = str_split($val, 11);
}
}
foreach($fca["hourly"] as $k => $fc){
$setval = trim($valA[$k]);
// ############ Aufbereitung der Daten je nach Daten aus XML ################
if($idstr == "SunD1") $setval = $setval / 60;
if($idstr == "RRad1") $setval = floatval($setval);
if($idstr == "T5cm") $setval = $setval - 273;
$fca["hourly"][$k][$idtxt] = $setval;
}
}
// Variable anlgen by Name
function CreateVariableByName($id, $name, $type, $profile = "")
{
# type: 0=boolean, 1 = integer, 2 = float, 3 = string;
global $_IPS;
$vid = @IPS_GetVariableIDByName($name, $id);
if($vid === false)
{
$vid = IPS_CreateVariable($type);
IPS_SetParent($vid, $id);
IPS_SetName($vid, $name);
IPS_SetInfo($vid, "this variable was created by script #".$_IPS['SELF']);
if($profile !== "") { IPS_SetVariableCustomProfile($vid, $profile); }
}
return $vid;
}
/* Sonnenstandsberechnung */
function calcSun($ts, $dLongitude, $dLatitude){
// Correction Timezone
$ts = $ts - 2*3600;
$iYear = date("Y", $ts);
$iMonth = date("m", $ts);
$iDay = date("d", $ts);
$dHours = date("H", $ts);
$dMinutes = date("i", $ts);
$dSeconds = date("s", $ts);
$pi = 3.14159265358979323846;
$twopi = (2*$pi);
$rad = ($pi/180);
$dEarthMeanRadius = 6371.01; // In km
$dAstronomicalUnit = 149597890; // In km
// Calculate difference in days between the current Julian Day
// and JD 2451545.0, which is noon 1 January 2000 Universal Time
// Calculate time of the day in UT decimal hours
$dDecimalHours = floatval($dHours) + (floatval($dMinutes) + floatval($dSeconds) / 60.0 ) / 60.0;
// Calculate current Julian Day
$iYfrom2000 = $iYear;//expects now as YY ;
$iA= (14 - ($iMonth)) / 12;
$iM= ($iMonth) + 12 * $iA -3;
$liAux3=(153 * $iM + 2)/5;
$liAux4= 365 * ($iYfrom2000 - $iA);
$liAux5= ( $iYfrom2000 - $iA)/4;
$dElapsedJulianDays= floatval(($iDay + $liAux3 + $liAux4 + $liAux5 + 59)+ -0.5 + $dDecimalHours/24.0);
// Calculate ecliptic coordinates (ecliptic longitude and obliquity of the
// ecliptic in radians but without limiting the angle to be less than 2*Pi
// (i.e., the result may be greater than 2*Pi)
$dOmega= 2.1429 - 0.0010394594 * $dElapsedJulianDays;
$dMeanLongitude = 4.8950630 + 0.017202791698 * $dElapsedJulianDays; // Radians
$dMeanAnomaly = 6.2400600 + 0.0172019699 * $dElapsedJulianDays;
$dEclipticLongitude = $dMeanLongitude + 0.03341607 * sin( $dMeanAnomaly ) + 0.00034894 * sin( 2 * $dMeanAnomaly ) -0.0001134 -0.0000203 * sin($dOmega);
$dEclipticObliquity = 0.4090928 - 6.2140e-9 * $dElapsedJulianDays +0.0000396 * cos($dOmega);
// Calculate celestial coordinates ( right ascension and declination ) in radians
// but without limiting the angle to be less than 2*Pi (i.e., the result may be
// greater than 2*Pi)
$dSin_EclipticLongitude = sin( $dEclipticLongitude );
$dY1 = cos( $dEclipticObliquity ) * $dSin_EclipticLongitude;
$dX1 = cos( $dEclipticLongitude );
$dRightAscension = atan2( $dY1,$dX1 );
if( $dRightAscension < 0.0 ) $dRightAscension = $dRightAscension + $twopi;
$dDeclination = asin( sin( $dEclipticObliquity )* $dSin_EclipticLongitude );
// Calculate local coordinates ( azimuth and zenith angle ) in degrees
$dGreenwichMeanSiderealTime = 6.6974243242 + 0.0657098283 * $dElapsedJulianDays + $dDecimalHours;
$dLocalMeanSiderealTime = ($dGreenwichMeanSiderealTime*15 + $dLongitude)* $rad;
$dHourAngle = $dLocalMeanSiderealTime - $dRightAscension;
$dLatitudeInRadians = $dLatitude * $rad;
$dCos_Latitude = cos( $dLatitudeInRadians );
$dSin_Latitude = sin( $dLatitudeInRadians );
$dCos_HourAngle= cos( $dHourAngle );
$dZenithAngle = (acos( $dCos_Latitude * $dCos_HourAngle * cos($dDeclination) + sin( $dDeclination )* $dSin_Latitude));
$dY = -sin( $dHourAngle );
$dX = tan( $dDeclination )* $dCos_Latitude - $dSin_Latitude * $dCos_HourAngle;
$dAzimuth = atan2( $dY, $dX );
if ( $dAzimuth < 0.0 )
$dAzimuth = $dAzimuth + $twopi;
$dAzimuth = $dAzimuth / $rad;
// Parallax Correction
$dParallax = ($dEarthMeanRadius / $dAstronomicalUnit) * sin( $dZenithAngle);
$dZenithAngle = ($dZenithAngle + $dParallax) / $rad;
$dElevation = 90 - $dZenithAngle;
return Array("azimuth" => $dAzimuth, "elevation" => $dElevation);
}
/*
* @author E. Akerboom
* @author {@link http://www.tremani.nl/ Tremani}, {@link http://maps.google.com/maps?f=q&hl=en&q=delft%2C+the+netherlands&ie=UTF8&t=k&om=1&ll=53.014783%2C4.921875&spn=36.882665%2C110.566406&z=4 Delft}, The Netherlands
* @since feb 2007
* @version 1.1
* @license http://opensource.org/licenses/bsd-license.php BSD License
*/
class NeuralNetwork {
protected $nodeCount = array ();
protected $nodeValue = array ();
protected $nodeThreshold = array ();
protected $edgeWeight = array ();
protected $learningRate = array (0.1);
protected $layerCount = 0;
protected $previousWeightCorrection = array ();
protected $momentum = 0.8;
protected $isVerbose = true;
protected $weightsInitialized = false;
public $trainInputs = array ();
public $trainOutput = array ();
public $trainDataID = array ();
public $controlInputs = array ();
public $controlOutput = array ();
public $controlDataID = array ();
protected $epoch;
protected $errorTrainingset;
protected $errorControlset;
protected $success;
/**
* Creates a neural network.
*
* Example:
* <code>
* // create a network with 4 input nodes, 10 hidden nodes, and 4 output nodes
* $n = new NeuralNetwork(4, 10, 4);
*
* // create a network with 4 input nodes, 1 hidden layer with 10 nodes,
* // another hidden layer with 10 nodes, and 4 output nodes
* $n = new NeuralNetwork(4, 10, 10, 4);
*
* // alternative syntax
* $n = new NeuralNetwork(array(4, 10, 10, 4));
* </code>
*
* @param array $nodeCount The number of nodes in the consecutive layers.
*/
public function __construct($nodeCount) {
if (!is_array($nodeCount)) {
$nodeCount = func_get_args();
}
$this->nodeCount = $nodeCount;
// store the number of layers
$this->layerCount = count($this->nodeCount);
}
/**
* Exports the neural network
*
* @returns array
*/
public function export()
{
return array(
'layerCount' => $this->layerCount,
'nodeCount' => $this->nodeCount,
'edgeWeight' => $this->edgeWeight,
'nodeThreshold' => $this->nodeThreshold,
'learningRate' => $this->learningRate,
'momentum' => $this->momentum,
'isVerbose' => $this->isVerbose,
'weightsInitialized' => $this->weightsInitialized,
);
}
/**
* Import a neural network
* @param array $nn_array An array of the neural network parameters
*/
public function import($nn_array)
{
foreach ($nn_array as $key => $value)
{
$this->$key = $value;
}
return $this;
}
/**
* Sets the learning rate between the different layers.
*
* @param array $learningRate An array containing the learning rates [range 0.0 - 1.0].
* The size of this array is 'layerCount - 1'. You might also provide a single number. If that is
* the case, then this will be the learning rate for the whole network.
*/
public function setLearningRate($learningRate) {
if (!is_array($learningRate)) {
$learningRate = func_get_args();
}
$this->learningRate = $learningRate;
}
/**
* Gets the learning rate for a specific layer
*
* @param int $layer The layer to obtain the learning rate for
* @return float The learning rate for that layer
*/
public function getLearningRate($layer) {
if (array_key_exists($layer, $this->learningRate)) {
return $this->learningRate[$layer];
}
return $this->learningRate[0];
}
/**
* Sets the 'momentum' for the learning algorithm. The momentum should
* accelerate the learning process and help avoid local minima.
*
* @param float $momentum The momentum. Must be between 0.0 and 1.0; Usually between 0.5 and 0.9
*/
public function setMomentum($momentum) {
$this->momentum = $momentum;
}
/**
* Gets the momentum.
*
* @return float The momentum
*/
public function getMomentum() {
return $this->momentum;
}
/**
* Calculate the output of the neural network for a given input vector
*
* @param array $input The vector to calculate
* @return mixed The output of the network
*/
public function calculate($input) {
// put the input vector on the input nodes
foreach ($input as $index => $value) {
$this->nodeValue[0][$index] = $value;
}
// iterate the hidden layers
for ($layer = 1; $layer < $this->layerCount; $layer ++) {
$prev_layer = $layer -1;
// iterate each node in this layer
for ($node = 0; $node < ($this->nodeCount[$layer]); $node ++) {
$node_value = 0.0;
// each node in the previous layer has a connection to this node
// on basis of this, calculate this node's value
for ($prev_node = 0; $prev_node < ($this->nodeCount[$prev_layer]); $prev_node ++) {
$inputnode_value = $this->nodeValue[$prev_layer][$prev_node];
$edge_weight = $this->edgeWeight[$prev_layer][$prev_node][$node];
$node_value = $node_value + ($inputnode_value * $edge_weight);
}
// apply the threshold
$node_value = $node_value - $this->nodeThreshold[$layer][$node];
// apply the activation function
$node_value = $this->activation($node_value);
// remember the outcome
$this->nodeValue[$layer][$node] = $node_value;
}
}
// return the values of the last layer (the output layer)
return $this->nodeValue[$this->layerCount - 1];
}
/**
* Implements the standard (default) activation function for backpropagation networks,
* the 'tanh' activation function.
*
* @param float $value The preliminary output to apply this function to
* @return float The final output of the node
*/
protected function activation($value) {
return tanh($value);
// return (1.0 / (1.0 + exp(- $value)));
}
/**
* Implements the derivative of the activation function. By default, this is the
* inverse of the 'tanh' activation function: 1.0 - tanh($value)*tanh($value);
*
* @param float $value 'X'
* @return $float
*/
protected function derivativeActivation($value) {
$tanh = tanh($value);
return 1.0 - $tanh * $tanh;
//return $value * (1.0 - $value);
}
/**
* Add a test vector and its output
*
* @param array $input An input vector
* @param array $output The corresponding output
* @param int $id (optional) An identifier for this piece of data
*/
public function addTestData($input, $output, $id = null) {
$index = count($this->trainInputs);
foreach ($input as $node => $value) {
$this->trainInputs[$index][$node] = $value;
}
foreach ($output as $node => $value) {
$this->trainOutput[$index][$node] = $value;
}
$this->trainDataID[$index] = $id;
}
/**
* Returns the identifiers of the data used to train the network (if available)
*
* @return array An array of identifiers
*/
public function getTestDataIDs() {
return $this->trainDataID;
}
/**
* Add a set of control data to the network.
*
* This set of data is used to prevent 'overlearning' of the network. The
* network will stop training if the results obtained for the control data
* are worsening.
*
* The data added as control data is not used for training.
*
* @param array $input An input vector
* @param array $output The corresponding output
* @param int $id (optional) An identifier for this piece of data
*/
public function addControlData($input, $output, $id = null) {
$index = count($this->controlInputs);
foreach ($input as $node => $value) {
$this->controlInputs[$index][$node] = $value;
}
foreach ($output as $node => $value) {
$this->controlOutput[$index][$node] = $value;
}
$this->controlDataID[$index] = $id;
}
/**
* Returns the identifiers of the control data used during the training
* of the network (if available)
*
* @return array An array of identifiers
*/
public function getControlDataIDs() {
return $this->controlDataID;
}
/**
* Shows the current weights and thresholds
*
* @param boolean $force Force the output, even if the network is {@link setVerbose() not verbose}.
*/
public function showWeights($force = false) {
if ($this->isVerbose() || $force) {
echo "<hr>";
echo "<br />Weights: <pre>".print_r($this->edgeWeight, true)."</pre>";
echo "<br />Thresholds: <pre>".print_r($this->nodeThreshold, true)."</pre>";
}
}
/**
* Determines if the neural network displays status and error messages. By default, it does.
*
* @param boolean $isVerbose 'true' if you want to display status and error messages, 'false' if you don't
*/
public function setVerbose($isVerbose) {
$this->isVerbose = $isVerbose;
}
/**
* Returns whether or not the network displays status and error messages.
*
* @return boolean 'true' if status and error messages are displayed, 'false' otherwise
*/
public function isVerbose() {
return $this->isVerbose;
}
/**
* Loads a neural network from a file saved by the 'save()' function. Clears
* the training and control data added so far.
*
* @param string $filename The filename to load the network from
* @return boolean 'true' on success, 'false' otherwise
*/
public function load($filename) {
if (file_exists($filename)) {
$data = parse_ini_file($filename);
if (array_key_exists("edges", $data) && array_key_exists("thresholds", $data)) {
// make sure all standard preparations performed
$this->initWeights();
// load data from file
$this->edgeWeight = unserialize($data['edges']);
$this->nodeThreshold = unserialize($data['thresholds']);
$this->weightsInitialized = true;
// load IDs of training and control set
if (array_key_exists("training_data", $data) && array_key_exists("control_data", $data)) {
// load the IDs
$this->trainDataID = unserialize($data['training_data']);
$this->controlDataID = unserialize($data['control_data']);
// if we do not reset the training and control data here, then we end up
// with a bunch of IDs that do not refer to the actual data we're training
// the network with.
$this->controlInputs = array ();
$this->controlOutput = array ();
$this->trainInputs = array ();
$this->trainOutput = array ();
}
return true;
}
}
return false;
}
/**
* Saves a neural network to a file
*
* @param string $filename The filename to save the neural network to
* @return boolean 'true' on success, 'false' otherwise
*/
public function save($filename) {
$f = fopen($filename, "w");
if ($f) {
fwrite($f, "[weights]");
fwrite($f, "\r\nedges = \"".serialize($this->edgeWeight)."\"");
fwrite($f, "\r\nthresholds = \"".serialize($this->nodeThreshold)."\"");
fwrite($f, "\r\n");
fwrite($f, "[identifiers]");
fwrite($f, "\r\ntraining_data = \"".serialize($this->trainDataID)."\"");
fwrite($f, "\r\ncontrol_data = \"".serialize($this->controlDataID)."\"");
fclose($f);
return true;
}
return false;
}
/**
* Resets the state of the neural network, so it is ready for a new
* round of training.
*/
public function clear() {
$this->initWeights();
}
/**
* Start the training process
*
* @param int $maxEpochs The maximum number of epochs
* @param float $maxError The maximum squared error in the training data
* @return bool 'true' if the training was successful, 'false' otherwise
*/
public function train($maxEpochs = 500, $maxError = 0.01) {
if (!$this->weightsInitialized) {
$this->initWeights();
}
if ($this->isVerbose()) {
echo "<table>";
echo "<tr><th>#</th><th>error(trainingdata)</th><th>error(controldata)</th><th>slope(error(controldata))</th></tr>";
}
$epoch = 0;
$errorControlSet = array ();
$avgErrorControlSet = array ();
$sample_count = 10;
do {
// echo "<tr><td colspan=10><b>epoch $epoch</b></td></tr>";
for ($i = 0; $i < count($this->trainInputs); $i ++) {
// select a training pattern at random
$index = mt_rand(0, count($this->trainInputs) - 1);
// determine the input, and the desired output
$input = $this->trainInputs[$index];
$desired_output = $this->trainOutput[$index];
// calculate the actual output
$output = $this->calculate($input);
// echo "<tr><td></td><td>Training set $i</td><td>input = (" . implode(", ", $input) . ")</td>";
// echo "<td>desired = (" . implode(", ", $desired_output) . ")</td>";
// echo "<td>output = (" . implode(", ", $output) .")</td></tr>";
// change network weights
$this->backpropagate($output, $desired_output);
}
// buy some time
#set_time_limit(300);
//display the overall network error after each epoch
$squaredError = $this->squaredErrorEpoch();
if ($epoch % 2 == 0) {
$squaredErrorControlSet = $this->squaredErrorControlSet();
$errorControlSet[] = $squaredErrorControlSet;
if (count($errorControlSet) > $sample_count) {
$avgErrorControlSet[] = array_sum(array_slice($errorControlSet, -$sample_count)) / $sample_count;
}
list ($slope, $offset) = $this->fitLine($avgErrorControlSet);
$controlset_msg = $squaredErrorControlSet;
} else {
$controlset_msg = "";
}
if ($this->isVerbose()) {
echo "<tr><td><b>$epoch</b></td><td>$squaredError</td><td>$controlset_msg";
echo "<script type='text/javascript'>window.scrollBy(0,100);</script>";
echo "</td><td>$slope</td></tr>";
echo "</td></tr>";
flush();
ob_flush();
}
// conditions for a 'successful' stop:
// 1. the squared error is now lower than the provided maximum error
$stop_1 = $squaredError <= $maxError || $squaredErrorControlSet <= $maxError;
// conditions for an 'unsuccessful' stop
// 1. the maximum number of epochs has been reached
$stop_2 = $epoch ++ > $maxEpochs;
// 2. the network's performance on the control data is getting worse
$stop_3 = $slope > 0;
} while (!$stop_1 && !$stop_2 && !$stop_3);
$this->setEpoch($epoch);
$this->setErrorTrainingSet($squaredError);
$this->setErrorControlSet($squaredErrorControlSet);
$this->setTrainingSuccessful($stop_1);
if ($this->isVerbose()) {
echo "</table>";
}
return $stop_1;
}
/**
* After training, this function is used to store the number of epochs the network
* needed for training the network. An epoch is defined as the number of times
* the complete trainingset is used for training.
*
* @param int $epoch
*/
private function setEpoch($epoch) {
$this->epoch = $epoch;
}
/**
* Gets the number of epochs the network needed for training.
*
* @return int The number of epochs.
*/
public function getEpoch() {
return $this->epoch;
}
/**
* After training, this function is used to store the squared error between the
* desired output and the obtained output of the training data.
*
* @param float $error The squared error of the training data
*/
private function setErrorTrainingSet($error) {
$this->errorTrainingset = $error;
}
/**
* Gets the squared error between the desired output and the obtained output of
* the training data.
*
* @return float The squared error of the training data
*/
public function getErrorTrainingSet() {
return $this->errorTrainingset;
}
/**
* After training, this function is used to store the squared error between the
* desired output and the obtained output of the control data.
*
* @param float $error The squared error of the control data
*/
private function setErrorControlSet($error) {
$this->errorControlset = $error;
}
/**
* Gets the squared error between the desired output and the obtained output of
* the control data.
*
* @return float The squared error of the control data
*/
public function getErrorControlSet() {
return $this->errorControlset;
}
/**
* After training, this function is used to store whether or not the training
* was successful.
*
* @param bool $success 'true' if the training was successful, 'false' otherwise
*/
private function setTrainingSuccessful($success) {
$this->success = $success;
}
/**
* Determines if the training was successful.
*
* @return bool 'true' if the training was successful, 'false' otherwise
*/
public function getTrainingSuccessful() {
return $this->success;
}
/**
* Finds the least square fitting line for the given data.
*
* This function is used to determine if the network is overtraining itself. If
* the line through the controlset's most recent squared errors is going 'up',
* then it's time to stop training.
*
* @param array $data The points to fit a line to. The keys of this array represent
* the 'x'-value of the point, the corresponding value is the
* 'y'-value of the point.
* @return array An array containing, respectively, the slope and the offset of the fitted line.
*/
private function fitLine($data) {
// based on
// http://mathworld.wolfram.com/LeastSquaresFitting.html
$n = count($data);
if ($n > 1) {
$sum_y = 0;
$sum_x = 0;
$sum_x2 = 0;
$sum_xy = 0;
foreach ($data as $x => $y) {
$sum_x += $x;
$sum_y += $y;
$sum_x2 += $x * $x;
$sum_xy += $x * $y;
}
// implementation of formula (12)
$offset = ($sum_y * $sum_x2 - $sum_x * $sum_xy) / ($n * $sum_x2 - $sum_x * $sum_x);
// implementation of formula (13)
$slope = ($n * $sum_xy - $sum_x * $sum_y) / ($n * $sum_x2 - $sum_x * $sum_x);
return array ($slope, $offset);
} else {
return array (0.0, 0.0);
}
}
/**
* Gets a random weight between [-0.25 .. 0.25]. Used to initialize the network.
*
* @return float A random weight
*/
private function getRandomWeight($layer) {
return ((mt_rand(0, 1000) / 1000) - 0.5) / 2;
}
/**
* Randomise the weights in the neural network
*/
private function initWeights() {
// assign a random value to each edge between the layers, and randomise each threshold
//
// 1. start at layer '1' (so skip the input layer)
for ($layer = 1; $layer < $this->layerCount; $layer ++) {
$prev_layer = $layer -1;
// 2. in this layer, walk each node
for ($node = 0; $node < $this->nodeCount[$layer]; $node ++) {
// 3. randomise this node's threshold
$this->nodeThreshold[$layer][$node] = $this->getRandomWeight($layer);
// 4. this node is connected to each node of the previous layer
for ($prev_index = 0; $prev_index < $this->nodeCount[$prev_layer]; $prev_index ++) {
// 5. this is the 'edge' that needs to be reset / initialised
$this->edgeWeight[$prev_layer][$prev_index][$node] = $this->getRandomWeight($prev_layer);
// 6. initialize the 'previous weightcorrection' at 0.0
$this->previousWeightCorrection[$prev_layer][$prev_index] = 0.0;
}
}
}
}
/**
* Performs the backpropagation algorithm. This changes the weights and thresholds of the network.
*
* @param array $output The output obtained by the network
* @param array $desired_output The desired output
*/
private function backpropagate($output, $desired_output) {
$errorgradient = array ();
$outputlayer = $this->layerCount - 1;
$momentum = $this->getMomentum();
// Propagate the difference between output and desired output through the layers.
for ($layer = $this->layerCount - 1; $layer > 0; $layer --) {
for ($node = 0; $node < $this->nodeCount[$layer]; $node ++) {
// step 1: determine errorgradient
if ($layer == $outputlayer) {
// for the output layer:
// 1a. calculate error between desired output and actual output
$error = $desired_output[$node] - $output[$node];
// 1b. calculate errorgradient
$errorgradient[$layer][$node] = $this->derivativeActivation($output[$node]) * $error;
} else {
// for hidden layers:
// 1a. sum the product of edgeWeight and errorgradient of the 'next' layer
$next_layer = $layer +1;
$productsum = 0;
for ($next_index = 0; $next_index < ($this->nodeCount[$next_layer]); $next_index ++) {
$_errorgradient = $errorgradient[$next_layer][$next_index];
$_edgeWeight = $this->edgeWeight[$layer][$node][$next_index];
$productsum = $productsum + $_errorgradient * $_edgeWeight;
}
// 1b. calculate errorgradient
$nodeValue = $this->nodeValue[$layer][$node];
$errorgradient[$layer][$node] = $this->derivativeActivation($nodeValue) * $productsum;
}
// step 2: use the errorgradient to determine a weight correction for each node
$prev_layer = $layer -1;
$learning_rate = $this->getlearningRate($prev_layer);
for ($prev_index = 0; $prev_index < ($this->nodeCount[$prev_layer]); $prev_index ++) {
// 2a. obtain nodeValue, edgeWeight and learning rate
$nodeValue = $this->nodeValue[$prev_layer][$prev_index];
$edgeWeight = $this->edgeWeight[$prev_layer][$prev_index][$node];
// 2b. calculate weight correction
$weight_correction = $learning_rate * $nodeValue * $errorgradient[$layer][$node];
// 2c. retrieve previous weight correction
$prev_weightcorrection = @$this->previousWeightCorrection[$layer][$node];
// 2d. combine those ('momentum learning') to a new weight
$new_weight = $edgeWeight + $weight_correction + $momentum * $prev_weightcorrection;
// 2e. assign the new weight to this edge
$this->edgeWeight[$prev_layer][$prev_index][$node] = $new_weight;
// 2f. remember this weightcorrection
$this->previousWeightCorrection[$layer][$node] = $weight_correction;
}
// step 3: use the errorgradient to determine threshold correction
$threshold_correction = $learning_rate * -1 * $errorgradient[$layer][$node];
$new_threshold = $this->nodeThreshold[$layer][$node] + $threshold_correction;
$this->nodeThreshold[$layer][$node] = $new_threshold;
}
}
}
/**
* Calculate the root-mean-squared error of the output, given the
* trainingdata.
*
* @return float The root-mean-squared error of the output
*/
private function squaredErrorEpoch() {
$RMSerror = 0.0;
for ($i = 0; $i < count($this->trainInputs); $i ++) {
$RMSerror += $this->squaredError($this->trainInputs[$i], $this->trainOutput[$i]);
}
$RMSerror = $RMSerror / count($this->trainInputs);
return sqrt($RMSerror);
}
/**
* Calculate the root-mean-squared error of the output, given the
* controldata.
*
* @return float The root-mean-squared error of the output
*/
private function squaredErrorControlSet() {
if (count($this->controlInputs) == 0) {
return 1.0;
}
$RMSerror = 0.0;
for ($i = 0; $i < count($this->controlInputs); $i ++) {
$RMSerror += $this->squaredError($this->controlInputs[$i], $this->controlOutput[$i]);
}
$RMSerror = $RMSerror / count($this->controlInputs);
return sqrt($RMSerror);
}
/**
* Calculate the root-mean-squared error of the output, given the
* desired output.
*
* @param array $input The input to test
* @param array $desired_output The desired output
* @return float The root-mean-squared error of the output compared to the desired output
*/
private function squaredError($input, $desired_output) {
$output = $this->calculate($input);
$RMSerror = 0.0;
foreach ($output as $node => $value) {
//calculate the error
$error = $output[$node] - $desired_output[$node];
$RMSerror = $RMSerror + ($error * $error);
}
return $RMSerror;
}
}
</code>
Forcast Ergebnis nach 3 Tagen Training (etwas schlechtes Wetter derzeit):