architecture-building-systems / CityEnergyAnalyst

The City Energy Analyst (CEA)
https://www.cityenergyanalyst.com/
MIT License
195 stars 65 forks source link

Demand calculation for HQ Masterplan scenario not working anymore #341

Closed martin-mosteiro closed 8 years ago

martin-mosteiro commented 8 years ago

After changing the occupancies and deleting the Shape_Leng and Shape_Area fields in the input files to be consistent with the baseline case, the demand calculation for the masterplan scenario is not working anymore...

Executing: DemandTool "C:\Users\User\Documents\GitHub\cea-reference-case\reference-case (HQ)\masterplan" Zurich
Start Time: Thu Sep 15 17:16:34 2016
Running script DemandTool...
Using python: C:\Python27\ArcGIS10.4\python.exe
Path to demand script: C:\Users\User\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py
Adding path to PYTHONPATH: C:\Users\User\Documents\GitHub\CEAforArcGIS
Running demand calculation for scenario C:\Users\User\Documents\GitHub\cea-reference-case\reference-case (HQ)\masterplan
Running demand calculation with weather file C:\Users\User\Documents\GitHub\CEAforArcGIS\cea\databases\CH\Weather\Zurich.epw
reading input files
done
calculating thermal properties
done
creating windows
done
Using 2 CPU's
Building No. 1 completed out of 242
Building No. 2 completed out of 242
Building No. 3 completed out of 242
Building No. 4 completed out of 242
Building No. 5 completed out of 242
Building No. 6 completed out of 242
Building No. 7 completed out of 242
Building No. 8 completed out of 242
Building No. 9 completed out of 242
Building No. 10 completed out of 242
Building No. 11 completed out of 242
Building No. 12 completed out of 242
Building No. 13 completed out of 242
Building No. 14 completed out of 242
Building No. 15 completed out of 242
Building No. 16 completed out of 242
Traceback (most recent call last):
 File "C:\Users\User\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 182, in <module>
   run_as_script(scenario_path=args.scenario, weather_path=args.weather)
 File "C:\Users\User\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 173, in run_as_script
   demand_calculation(locator=locator, weather_path=weather_path, gv=gv)
 File "C:\Users\User\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 105, in demand_calculation
   weather_data)
 File "C:\Users\User\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 153, in thermal_loads_all_buildings_multiprocessing
   job.get(240)
 File "C:\Python27\ArcGIS10.4\lib\multiprocessing\pool.py", line 567, in get
   raise self._value
RuntimeError: Failed to converge after 100 iterations, value is (nan+nan*j)

Completed script DemandTool...
Succeeded at Thu Sep 15 17:17:11 2016 (Elapsed Time: 37.07 seconds)
daren-thomas commented 8 years ago

I was able to reproduce the bug... I will look into fixing it on monday...

daren-thomas commented 8 years ago

BTW: setting multiprocessing to False makes the error message clearer:

Traceback (most recent call last):
  File "C:\Users\darthoma\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 182, in <module>
    run_as_script(scenario_path=args.scenario, weather_path=args.weather)
  File "C:\Users\darthoma\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 173, in run_as_script
    demand_calculation(locator=locator, weather_path=weather_path, gv=gv)
  File "C:\Users\darthoma\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 108, in demand_calculation
    weather_data)
  File "C:\Users\darthoma\Documents\GitHub\CEAforArcGIS\cea\demand\demand_main.py", line 138, in thermal_loads_all_buildings
    building, bpr, weather_data, usage_schedules, date, gv, locator)
  File "C:\Users\darthoma\Documents\GitHub\CEAforArcGIS\cea\demand\thermal_loads.py", line 257, in calc_thermal_loads
    tsd['ta_hs_set'])
  File "C:\Users\darthoma\Documents\GitHub\CEAforArcGIS\cea\demand\sensible_loads.py", line 351, in calc_temperatures_emission_systems
    Ta_sup_0, Ta_re_0, gv.Cpa)
  File "C:\Python27\ArcGIS10.4\lib\site-packages\numpy\lib\function_base.py", line 1700, in __call__
    return self._vectorize_call(func=func, args=vargs)
  File "C:\Python27\ArcGIS10.4\lib\site-packages\numpy\lib\function_base.py", line 1769, in _vectorize_call
    outputs = ufunc(*inputs)
  File "C:\Users\darthoma\Documents\GitHub\CEAforArcGIS\cea\technologies\heating_coils.py", line 86, in calc_cooling_coil
    result = sopt.newton(fh, trc0, maxiter=100, tol=0.01) - 273
  File "C:\Python27\ArcGIS10.4\lib\site-packages\scipy\optimize\zeros.py", line 161, in newton
    raise RuntimeError(msg)
RuntimeError: Failed to converge after 100 iterations, value is (nan+nan*j)
daren-thomas commented 8 years ago

The building in question is: MP09

daren-thomas commented 8 years ago

Here is some more information: The error happens in the newton solver for the cooling coil. At the time of error, inside the function fh, the following variables are set:

'Qcsf_0': -2160935.2513834536, 
'mCw0': 270116.9064229317, 
'LMRT': (nan+nan*j), 
'k2': 0.5015153719375163, 
'x': 288L, 
'LMRT0': (-1.0775220302930815-1.9515353596348863j), 
'TD2': 0j, 
'TD1': (-0.50151537193750073+0j), 
'Eq': (nan+nan*j), 
'tc': (288+0j)}

The main problem seems to be the line LMRT = (TD2 - TD1) / scipy.log(TD2 / TD1), as the value of TD2/TD1 results in (-0-0j) and well, the log for that is not a happy number: (-inf-3.1415926535897931j). And a runtime warning (scipy.log(TD2 / TD1) C:\Users\darthoma\Anaconda2\envs\esri104\lib\site-packages\numpy\lib\scimath.py:262: RuntimeWarning: divide by zero encountered in log)

@JIMENOFONSECA, @gabriel-happle: We can check for this runtime warning, no problem. But what should be the output in such a case? That is, how should we handle this problem? I think this is for domain experts to do...

daren-thomas commented 8 years ago

OK. I'm working on a fix that uses the bisect algorithm as a fallback when the newton fails.

daren-thomas commented 8 years ago

Check out PR #342 - it includes a fallback to bisect (slower but more reliable) if the newton fails for the cooling coil temperature calculation. The masterplan runs through, but with a ComplexWarning (the complex part of the calculation seems to be truncated, but I don't know if this is a problem at all...). It would be good to check the results using @gabriel-happle 's debug output method.

jimenofonseca commented 8 years ago

Hi, in my opinion there is no need of the bisct algorithm. This error has happened in the past. The issue is that the inputs of the function are just wrong. Example: negative temperatures, deltaT=0. Please take a look at the value of the inputs not at the newton function. There is nothing wrong with it...

Sorry for my late response for this issue. As discussed, we were away for some time...

On 19 Sep 2016, at 21:48, Daren Thomas notifications@github.com<mailto:notifications@github.com> wrote:

Check out PR #342https://github.com/architecture-building-systems/CEAforArcGIS/pull/342 - it includes a fallback to bisect (slower but more reliable) if the newton fails for the cooling coil temperature calculation. The masterplan runs through, but with a ComplexWarning (the complex part of the calculation seems to be truncated, but I don't know if this is a problem at all...). It would be good to check the results using @gabriel-happlehttps://github.com/gabriel-happle 's debug output method.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/architecture-building-systems/CEAforArcGIS/issues/341#issuecomment-247997847, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AIjrgnxNlGL3CJWJYvIQzKpsnSPkCbUcks5qrpKrgaJpZM4J-AKi.

daren-thomas commented 8 years ago

In this case, the documentation of the function should state restrictions on inputs, e.g. a valid range of inputs. Also, unit testing should test the boundaries of this function. I will create an issue to fix this. @JIMENOFONSECA Maybe we can work together on that one? This is probably also the case for the heating coil and tabs calculations.

Ideally, a function should check its inputs for valid ranges and throw an error if they are not met. The unit tests test to make sure this happens. I am currently working on installing a CI server (Jenkins is high on the list) on my old PC that will run these tests every night / on every commit to master.

daren-thomas commented 8 years ago

We are still left with the problem that the program creates an invalid state - the bisect algorithm patched that under the assumption that the newton algorithm does not converge for all cases - which is generally true but might not be in this special case.

@JIMENOFONSECA: I think we might need help debugging this problem as your domain knowledge would greatly increase our chances of finding the bug quickly.

daren-thomas commented 8 years ago

This is the list of local variables when the newton fails:

'Qcsf_0': -2160935.2513834536, 
'ec': (0.99496264113301858+0.071115583768434992j), 
'Ta_re_cs': 15.0, 
'Ta_sup_0': 281.2, 
'Qcsf': -135467.78079130786, 
'Ta_sup_cs': 15.0, 
'tare': 288.0, 
'k2': 0.5015153719375163, 
'UA0': (468545.07095939072-848597.28882868378j), 
'Tcs_sup_0': 7L, 
'tc': (288+0j), 
'tsc0': 280L, 
'TD10': -6.800000000000011, 
'ma_sup_0': 392.3009555500433, 
'mCw0': 270116.9064229317, 
'Cpa': 1.008, 
'trc0': 288L, 
'Ta_re_0': 281.2, 
'tasup': 288.0, 
'AUa': (32019.380751591805-57991.34678791475j), 
'NTUc': (2.6409463721774475-4.7831042737956126j), 
'ma_sup_cs': 12.02798323236849, 
'Tcs_re_0': 15L, 
'LMRT0': (-1.0775220302930815-1.9515353596348863j), 
'TD20': 1.1999999999999886

The inputs to the function are: Qcsf, Qcsf_0, Ta_sup_cs, Ta_re_cs, Tcs_sup_0, Tcs_re_0, ma_sup_cs, ma_sup_0, Ta_sup_0, Ta_re_0, Cpa

jimenofonseca commented 8 years ago

the problem is with Ta_re_cs and Ta_sup_cs being the same. does this happen for the time before? are we using the version of CEA of the master right?

daren-thomas commented 8 years ago

I moved to the branch of the pull merge, but it is up-to-date with master. Yes, I saw these are the same. I am now testing to see if there are other instances when they are the same.

How can they end up that way?

jimenofonseca commented 8 years ago

not sure... but the best way to see it is to check the value of the variables before the error. did you check the inputs for the building in the shapefiles? It can be that some of the properties are 0 and then the algorithm just crashes when it does not find a balance.

daren-thomas commented 8 years ago

Ta_re_cs and Ta_sup_cs are often equal! I just checked. In fact, when running the masterplan scenario, this state happens 10'860 times! And it does so for a lot of different buildings...

It only falls back to the bisect once.

jimenofonseca commented 8 years ago

mmm....ok i will need to see the data... coudl you upload the case study and tell me which building name?

daren-thomas commented 8 years ago

The case study is the hq masterplan from the repository. The building name is MP09.

martin-mosteiro commented 8 years ago

@JIMENOFONSECA: The problem seems to be the new buildings, not the old ones. out of these, the only values in the input files that are zero are height below ground (for some buildings which have none) and renovation years - which are zero since they have never been renovated. So unless the renovation year convention has changed (I can't seem to find it in the manual anymore), the values given should be correct.

martin-mosteiro commented 8 years ago

All other inputs (including those generated by the data helper) appear reasonable.

jimenofonseca commented 8 years ago

what is the occupancy of the building?

martin-mosteiro commented 8 years ago

Hospital

martin-mosteiro commented 8 years ago

The scenario works perfectly with Zug weather.

daren-thomas commented 8 years ago

I went ahead and plotted fh(x) for the offending case (see input data above):

image

Notice how there are multiple solutions? Also, note that python throws some warnings:

C:\Users\darthoma\Anaconda2\lib\site-packages\numpy\lib\scimath.py:262: RuntimeWarning: divide by zero encountered in log return nx.log(x) C:\Users\darthoma\Anaconda2\lib\site-packages\ipykernel\__main__.py:4: RuntimeWarning: invalid value encountered in cdouble_scalars

daren-thomas commented 8 years ago

Compare to a case with Ta_re_cs = Ta_sup_cs + 1:

image

daren-thomas commented 8 years ago

image

Note here that the starting value for x is the red line. A solver might run into trouble here. bisect finds this solution:

image

martin-mosteiro commented 8 years ago

Given that the inputs all seem to be consistent and the scenario works with other weather, I do think the multiple minimums that @daren-thomas are the most likely explanation. I would go forward with using the bisect algorithm for the MIBS course. I will continue testing the scenarios ( #340 ) using this fallback solution for now.

martin-mosteiro commented 8 years ago

The merge solves this issue. Since we already had similar issues with other weather files (Cairo and Sweden come to mind), maybe @amraladdin could test these cases again to see if this merge solves those previous issues as well.

However, the issue with the masterplan scenario has been solved, so it will now be closed.

jimenofonseca commented 8 years ago

Be aware you need to test the right weather file, so the right latitude and longitud. you have also to specify the cooling and heating season variable in global vars...

On 22 Sep 2016, at 3:55 pm, martin-mosteiro notifications@github.com<mailto:notifications@github.com> wrote:

The merge solves this issue. Since we already had similar issues with other weather files (Cairo and Sweden come to mind), maybe @amraladdinhttps://github.com/amraladdin could test these cases again to see if this merge solves those previous issues as well.

However, the issue with the masterplan scenario has been solved, so it will now be closed.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/architecture-building-systems/CEAforArcGIS/issues/341#issuecomment-248835972, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AIjrgnh6UwM3QWb1GlOj96PaLXNHKeU8ks5qsjRXgaJpZM4J-AKi.