Add a modular heat pump approach

[FWuellhorst]

Many researches need models of heat pumps and chillers, and there are various approaches with different levels of detail.
Currently, the IPBSA supports basically two approaches, leading to four models:

• WaterToWater models for scroll and reciproting compressor
• Carnot based appraoches, on with y and one with TCon as an input

The Carnot models are also available as a chiller option.
While the WaterToWater models explicitly model the refrigerant properties, many data driven / black-box approaches for the vapour compression cycle exists in literature. The Carnot model may be viewed as such a black-box approach.

To enable different black-box approaches for refrigerant cycles and thus heat pumps and chillers, we developed an approach in the AixLib and presented the results in the latest Modelica Conference (https://ecp.ep.liu.se/index.php/modelica/article/view/234).
The motivation for this model and the safety controls is stated in detail in the paper.
If there is an interest to integrate this approach into the core library, several steps would be necessary:

• Add the PartialVapourCompressionMachine and the extending Chiller and HeatPump model
• Add the SafetyControl models
• Add the validation data from the paper
• Add an example showing the interaction between heat pump and safety controls
• Revisit the implementation of the carnot models and check if both approaches can be merged. Currently, I only see issues regarding the TCon approach, as our approach takes the compressor frequency as an input, and the condenser outlet is an output of the model. However, we could think about some option to enable both control types, compressor speed and supply temperature. The last one would be especially relevant for certain MPC approaches.
• Check if the approach is comatible with the IBPSA standards, OpenModelica and can match the current approaches with regard to simulation speed

@mwetter : As we've already briefly disussed this in the last project 1 meeting. What is your opinion on this topic?

[mwetter]

I think this would be a good addition. Once it is added, we can see if it should indeed make the existing Carnot model obsolete, or whether both should be kept as the Carnot model is probably as simple as one can model a HP while taking the lift into account. I could see applications for early conceptual design/technology screening or MPC in which a simpler model that has not the discontinuities introduced by the safety control may still be applicable.

[FWuellhorst]

Thanks for the quick answer! I will start to work on this.
Regarding the discontinuities: The safety option is currently fully separated from the physical model to separate control and physics. Thus, the basic model has no safety control whatsoever.
But I think adding it first and checking if the Carnot-Approaches would be deprecated in future versions makes sense. Especially because the conversion scripts could not handle some change like this. At least I don't know of such commands. Removal would thus be non-backward-compatible.
Regarding the TCon approach: I just came across a use-case where I had to modify this in my approach. I think I will directly add an option for this when merging into IBPSA.

[FWuellhorst]

Additional tasks when integrating the model:

• Move the replaceable iceFac models into the physical model itself. This way, controls and physics are clearly separated
• Add TConOutSet option to enable control of the flow temperature. This will require a onOffSet as well.
• Automatic design: Change the data records to enable a design based on bivalence temperatures. This should be an advanced option, but an option nonetheless.
• Enable a replaceable inertia model. Currently, either the damper or no inertia is modeled. Using a constrainedby for a RealSISO interface, one could add for instance a delay to model dead bands.
• Consider finding a better name for the PartialPerformanceData. E.g.: ParialVapourCompressionBlackBoxModel,
• Move reduction of QEva_flow and QCon_flow into the innerCycle model to avoid having to implement this feature in all new PerformanceData approaches.

[FWuellhorst]

• Rename nSet to ySet
• Use hysteresis for onOffMea signal to avoid numerical issues

[FWuellhorst]

• Move Partial model to HeatPumps package
• Rename PerformanceData to BlackBox
• Align models with the default grid
• State the options for outputs only or outputs / inputs for the bus connectors
• Disable sensors and use inStream and volume ports directly
• Move regression tests from AixLib

[FWuellhorst]

• Move the scaling factor to the BlackBox data models
• Re-label the automatic design estimation -> Was only validated for industrial heat pumps?

[FWuellhorst]

@DaJansenGit: Could you check for which size-range your automatic design is valid? I would like to add a warning and maybe state more clearly, that it was done for large-scale heat pumps.

[DaJansenGit]

The data basis is ~ 20 data sets for water-to-water heat pumps in the range from 25 kW to 1 MW. This is already stated in the model information. However, I agree that an assert warning would help to avoid confusion.

[FWuellhorst]

@mwetter : One question regarding warnings during check and translation.
I added a method for an automatic design. This way, users can set nominal values as in the Carnot_ models, e.g. QCon_flow_nominal, TCon_nominal, TEva_nominal, ...
In the BlackBox-Data models, these values are used to calculate the nominal heat flow rate QConBlackBox_flow_nominal. For table-based data, I interpolate just once in translation. Using both values, the scalingFactor is calculated:

final parameter Real scalingFactor=QCon_flow_nominal/QConBlackBox_flow_nominal

In the AixLib version, this process is manual and very error-prone.
QConBlackBox_flow_nominal is thus a non-literal value. As the scalingFactor scales the nominal mass flow rates, m_flow_nominal becomes non-literal as well:

In the examples, I overwrite the method to calculate QConBlackBox_flow_nominal and set the literal value to avoid the warnings and thus fail of CI. However, this would not be the intended use of the model.
Any thoughts on what option is best? Is it important to mitigate the warnings besides the errors in the CI?

[FWuellhorst]

• Add minimal flow rate to security control

[FWuellhorst]

• Add validation package for data from Modelica Conference Paper

[mwetter]

@FWuellhorst : I would leave the warnings about the nominal attribute being ignored. It would be good to raise this problem to Dassault. We also have a few models in the Buildings library that exhibit exactly this issue. Best would be if Dymola can assign these nominal values once the model is initialized. Currently, they are written into dsmodel.c and hence the above warning needs to be issued.

[Mathadon]

Same in IDEAS. It would actually be even better if Dymola can assign the nominal values during translation. It should be perfectly possible to evaluate the Media functions during translation (unless their inputs depend on an initial equation). I think this actually works in OMC/JModelica but I'm not sure.

[FWuellhorst]

Ok, thanks for the feedback! Then I will continue with this new design approach.

@Mathadon : I checked and I don't get the error with the snippet below in OpenModelica:

model NonLiteralWarning
  parameter Real m_flow_nominal = if u then 1 else 0;
  Real m_flow(nominal=m_flow_nominal);
  parameter Boolean u annotation(Evaluate=false);
equation 
  m_flow = 0;
end NonLiteralWarning;

Dymola 2022x:

OpenModelica :

[FWuellhorst]

• Add a string identifier for data records and warn about diverging chiller and heat pump models

[FWuellhorst]

@mwetter @Mathadon :
Even though I had high hopes due to the successful check, the simulation in OpenModelica fails for the heat pump model as soon as the mass flow rate is non-literal in Dymola. Errors are raised that the mass flow rates are 0, and that m_flow_pos must be greater than 0.

Analyzing why this happens, I found the following: The scaling factor is calculated using an expression like this QUse/QUseBlackBoxData. If I add a start value to the parameter QUse, the simulation succeeds. However, the results are wrong. I propagate QUse from the uppermost level to the replaceable model, and OpenModelica just disregards this value and all other values given by the uppermost level. Maybe the replaceable model propagation does not work in OpenModelica. I followed this guide: https://www.claytex.com/blog/replaceable-model-declaration-propagation/

As I find this replaceable model definition in various locations in the MLS, I assume it's fine to use it (https://modelica.org/documents/MLS.pdf).

Here is a minimal example showing the issue (the heat pump approach uses this replaceable structure).

package TestModelPropagation

partial model PartialSubsystem
  parameter Real p;
  Real y;
end PartialSubsystem;

model QuadraticP
  extends PartialSubsystem;
  
  equation
  y = p * p;
end QuadraticP;

model UppermostSystem
  parameter Real p2;
  replaceable model MySubsystem = PartialSubsystem constrainedby PartialSubsystem(p=p2);
  MySubsystem subsys;
  Real y;
  
  equation
    y = subsys.y;
  
end UppermostSystem;

model Example

  UppermostSystem sys(
    p2=100,
    redeclare model MySubsystem = QuadraticP);
end Example;
  
end TestModelPropagation;

sys.y should be 100^2, but is zero.

I hope that OM compatibility is not a limiting factor when considering this model approach. I'm happy to raise this issue to the developers to fix this issue, or maybe we find a solution within my code that causes this error in the first place.

[Mathadon]

@FWuellhorst if there is a bug in OpenModelica than we can best report it there and hope that it will be fixed soon. So far their response has been rather quick so I have good hopes :)

[FWuellhorst]

The response was really fast, and the bug was already fixed in the latest version. The model works in OM 1.19.2 and simulates. At a first glance, the results at least look similar.

[mwetter]

@FWuellhorst Thanks for the minimum test model. I confirm that this works in OpenModelica 1.20.0~dev-250-gb17e1a0.
I am in the middle of updating OpenModelica for the CI tests (currently testing it for Buildings and then updating it for IBPSA as well) and afterwards this model should work too.

[FWuellhorst]

Thanks for checking and confirming. I will raise the compatibility issue in the AixLib as well :)

[Mathadon]

@mwetter fyi we have plans to use OM for CI testing in the future too.

[FWuellhorst]

Digging deeper into the differences between Dymola and OM for my model, I found an error in one safety control (algebraic loop) which shouldn't be possible in the first place. Basically, Dymola accepted ySet = ySet as a connection.
So it really seems that OM works as long as the models use proper Modelica.

One question regarding the CI pipelines: Are you using the code from buildingspy? I've seen that OM itself tests our libraries as well (https://libraries.openmodelica.org/branches/master/). Do you you their code for testing as well?

[mwetter]

@FWuellhorst We run our own tests with OpenModelica to catch errors before they show up in the CI testing of openmodelica.org.
See

modelica-ibpsa/.travis.yml

Line 24 in 2584136

- OMC_VERSION=ubuntu-2004-omc:1.19.0_dev-551-g1772c4b-1

for the specification of the OpenModelica image.

[FWuellhorst]

As discussed in #1628, the naming HeatPump and Chiller are not fitting with the usual approach in the IBPSA. The name should specify what the model is about.

I am currently in favour of ModularReversible, as it's the title of the paper. Adding BlackBox to the title is not fitting, as the basic model equations tend towards white-box. At the same time reversible could imply thermodynamically reversible, hence another Carnot approach. But I don't know a better term for a device which is able to both act as chiller and heat pump.

@DaJansenGit @LauraMaier : As users of the model in the AixLib, what do you think?

[FWuellhorst]

Based on the discussions with @mwetter @hcasperfu, the following new todos arise:

• Make the bus-connector non-visible or protected using a boolean option
• disable TSet if use_busConnectorsOnly is active
• rename mode to revert
• In the safety control blocks, change the use_ blocks to the upper level
• don't use mode apart from defrost
• All flows from left to right
• Check hysteresis in boundary map model
• Move automatic parameter estimation to a specific model for large scale water to water devices
• Change V to tau to follow IBPSA guidelines

[hcasperfu]

@FWuellhorst It was nice talking to you this morning in the meeting!
Just to clarify:

• All flows from left to right

This is (mainly) just for the control blocks.

I will also be happy to look at the inheritance structure of the code some time later to see if some of the code repetition can be reduced.

[FWuellhorst]

@mwetter : I used the following solution to disable the bus connector for Buildings and enable it in the AixLib:

  final parameter Boolean use_extBusCon=true
    "=false to disable external use of the bus connector"
    annotation(choices(checkBox=true), Dialog(group="Input Connectors"));

  Interfaces.VapourCompressionMachineControlBus sigBus if use_extBusCon
    "Bus with signal for device control" annotation (
          Placement(transformation(extent={{-120,-60},{-90,-26}}),
        iconTransformation(extent={{-108,-52},{-90,-26}})));
protected
  Interfaces.VapourCompressionMachineControlBus sigBusInt if not use_extBusCon
    "Internal sigBus to avoid external access";
equation
  // Internal bus connections --> No graphics required
  connect(mFlow_eva.m_flow, sigBusInt.m_flowEvaMea);
  // External bus connections
  connect(mFlow_eva.m_flow, sigBus.m_flowEvaMea) annotation(...);

Internally, sigBus is still used for the black box models.
Buildings could set the final parameter to false, AixLib to true.

One question remains regarding the default of this parameter in the IBPSA. The example would require to disable this option. Otherwise, you would have to exclude these examples upon merging into Buildings. Thus, the default false would make sense. At the same time, this would inhibit showing the other libraries how to use the connector for easy control.

[FWuellhorst]

• Rename to ModularReversible

• Add examples for replacing single points. Suggestions below:

• ModularReversible --> General modular model

• LargeScaleWaterToWater --> Model with automatic parameter estimation

• ReversibleAirToWaterEuropeanNorm2D --> Show 2D data and reversible operation

• ReversibleCarnotWithLosses --> Show usage of frost and heat losses to ambient

• AirToWaterWithFlowsheetAndFluid --> Show detailed white-box physics performance map

For each of these, an example will be introduced. With all 5, the usage should be clear enough. What do you think @hcasperfu @mwetter ?

[mwetter]

@FWuellhorst The list of examples look good.

Regarding the connector, I would rather not include examples in IBPSA that use the external signal bus, otherwise this may lead to a proliferation of models that uses this modeling concept. You could always add these examples in AixLib for the AixLib users. Would that work?
I am however not sure how to set the parameter final during the merge. This would need changes to BuildingsPy (at which point we could as well consider deleting what is always removed?).

[FWuellhorst]

Adding examples in the AixLib would work. I think I would just switch out the control blocks of the IBPSA examples.

Regarding the proliferation: It seemed to me that a lot of commercial libraries use expandable connectors, especially for control applications. Showing users that buses can be used aside from the weather information would, in my opinion, not be a bad thing.
But I see that no other models use this concept, and that it's also not compatible to CDL, or? I will set the default to false.

Regarding merge, maybe a solution like this would come in handy:

PartialVapourCompressionMachine: Define base components etc.
PartialVapourCompressionMachineWithBus (extends: PartialVapourCompressionMachine): Only connect the internal / external bus
ModularReversibleHeatPump and ModularReversibleChiller (all extend PartialVapourCompressionMachineWithBus): Define devices.

If all the relevant information is in PartialVapourCompressionMachineWithBus regarding the bus-connectors, the script in BuildingsPy could easily be adapted by a case just for that file. Whether you then change a boolean or remove the external bus is up to you.
What do you think?

[mwetter]

CDL does not allow expandable connectors, as the tooling for supporting them is quite complex.

I think BuildingsPy could be updated to delete certain lines during the merge. That probably leads to the cleanest end result.

[FWuellhorst]

Like this?

// Line to delete if you don't want to use the bus externally
protected 
  Interfaces.VapourCompressionMachineControlBus sigBus
    "Bus with signal for device control" annotation (
          Placement(transformation(extent={{-120,-60},{-90,-26}}),
        iconTransformation(extent={{-108,-52},{-90,-26}})));

  parameter Boolean use_busConnectorOnly=false
    "=true to use bus connector for model inputs (revSet, ySet, TSet, onOffSet).
    =false to use the bus connector for outputs only. 
    Only possible if no internal safety control is used"
    annotation(choices(checkBox=true), Dialog(group="Input Connectors", enable=not 
          use_safetyControl));

// Line to add if you want to use the bus
//protected
  parameter Real scaFac "Scaling-factor of vapour compression machine";
  parameter MediumCon.ThermodynamicState staCon_nominal=MediumCon.setState_pTX(
      T=MediumCon.T_default, p=MediumCon.p_default, X=MediumCon.X_default)
      "Nominal / default state of condenser medium";

  parameter MediumEva.ThermodynamicState staEva_nominal=MediumEva.setState_pTX(
      T=MediumEva.T_default, p=MediumEva.p_default, X=MediumEva.X_default)
      "Nominal / default state of evaporator medium";

equation

Libraries using the concept can remove the protected line and add it two lines below. Contrary, IBPSA defaults with no bus, and thus, Buildings as well.

[FWuellhorst]

As discussed in the IBPSA Working Group meeting, I reverted some features:

• SDF dependency
• TSet feature, tbd. in a future issue
• revSet is now hea and coo for heat pump and chiller, respectively.
• I moved the defrost safety block out of the safety control. It highly depends on the defrost model in use. Future issues should implement experimentally validated approaches, or at least literature based as by Zhu et al.
• I added the safety block also to the chiller. Chiller as well contain minimal runtime etc. It made no sense to just include this feature for heat pumps.
• I renamed Black-Box-Data and vapour compression cycle to refrigerant cycle. Vapour compression would be the more general term, also applying to power plants using coal etc.

[FWuellhorst]

Regarding the naming of not useful:

Current implementation

Heat exchangers are called con and eva, nominal values TCon_nominal etc., measured values TConOutMea etc.
Issues with that:

• If operation changes (cooling / heating), in reality, the condenser is used as an evaporator and vice versa
• Users have to constantly think about which side is what for reversible devices
• Operational Envelopes models require the user-side temperatures to check limits. If the temperatures are called TEva and TCon, it won't work on reverse mode. Currently, the connection switches the two temperatures and connects them to a TUse and TNotUse. The name NotUse is not really pretty.

The same issue holds for names like source/sink or load/source depending on if you look from the system or heat pump perspective.

Proposed Solution

@hcasperfu came up with ambient side AmbSid and user side UseSid. While not all use-cases will have human users, the 3-letter naming restriction actually helps in this case, as use is valid for all applications.
The tricky part is ambient, as it may be confused with the ambient air. For ground-source heat pumps, this would be confusing. Ambient side however, should be clear in all cases, as it is the ambient side of the heat pump / chiller. This also makes sense from a teaching perspective, as heat pumps are explained taking heat from the ambient, chillers exhausting the heat to the ambient.

Discussion

I could just use these names for the operational envelope. However, the concept could also be applied for con and eva. We would then have TAmbSid_nominal, TUseSid_nominal, TUseSidOutMea, etc. The heat flows would become QUseSid_flow, and the heat exchangers useSidHeaExc and ambSidHeaExc. The latter are not as straight forward as con and eva, but again, eva and con are not technically correct.

As the renaming would be some effort, let's decide before going into action. From my point of view, as the model is for system application, the naming could make the life of system engineers much simpler. At the same time, the existing models use con and eva as they are not reversible. What do you think @mwetter @hcasperfu ?

[mwetter]

We also had multiple iterations around that and could not come up with a good scheme. The issue of useful/not useful is that it won’t work for a heat recovery chiller, as both sides are useful.
For the operational envelop check, maybe we need two instances for this check, one for the normal operation and one for the reverse mode. I have not looked at the literature but I would think that the envelopes are different for the different modes?

Maybe for reversible heat pumps we need to get back and use a and b. That may also be a case where we need to use a different base class rather than trying to fit everything onto one base class and then having to use notations that are not intuitive to the user for a large class of applications.

[FWuellhorst]

For the operational envelope, I have two instances for main and second mode. I could make two operational envelope models, one for chillers and one for heat pumps. Having a common base, extending and just adding the input connectors with the correct names would avoid code repetition.

Regarding naming everything a, b or 1 and 2: This just shifts the naming issue as the user still has to check which side is the condenser / evaporator. Currently, the documentation clearly states that the heat exchanger names are for the main operation. If one models a heat pump, I think it's much clearer to have con and eva.
That being said, I would keep the names and just duplicate the operational envelope, ok?

[hcasperfu]

I agree with @FWuellhorst that naming things a and b is less descriptive and less helpful.
In my opinion the heat recovery chiller is less of an issue as it exceeds the intended scope of this model, even though the said chiller is perfectly within this model's functionality. What we want here is a clear naming scheme that accommodates the switchable heat pump model, which is what this model is advertised to do.

[mwetter]

Let's leave it at eva and con and be clear that this is for the main operation. I think that would be easiest for the user.

[FWuellhorst]

Just to note it here, a comparison between the modular models and the Carnot_y model in terms of model statistics and execution time.
I tested the execution time for 100 consecutive run with a pre-translated model using the experimentSetup settings.

	Carnot_y	ConstantQualityGrade	TableData2D
Original Model
Number of components	58	93	96
Variables	661	868	889
Constants	18 (18 scalars)	16 (16 scalars)	16 (16 scalars)
Parameters	311 (295 scalars)	423 (413 scalars)	454 (483 scalars)
Unknowns	332 (252 scalars)	429 (425 scalars)	419 (415 scalars)
Differentiated variables	4 scalars	4 scalars	4 scalars
Equations	272	420	414
Nontrivial	228	351	347
Translated Model
Constants	395 scalars	590 scalars	583 scalars
Free parameters	4 scalars	4 scalars	5 scalars
Parameter depending	26 scalars	26 scalars	72 scalars
Continuous time states	2 scalars	2 scalars	2 scalars
Time-varying variables	49 scalars	47 scalars	42 scalars
Alias variables	91 scalars	187 scalars	212 scalars
Number of mixed real/discrete systems of equations	0	0	0
Sizes of linear systems of equations	{ }	{ }	{ }
Sizes after manipulation of the linear systems	{ }	{ }	{ }
Sizes of nonlinear systems of equations	{10}	{8}	{ }
Sizes after manipulation of the nonlinear systems	{1}	{2}	{ }
Number of numerical Jacobians	0	0	0
Integration terminated successfully at T = 3600
Execution time for 100 simulation runs	42 seconds	42 seconds	41 seconds
CPU-time for initialization	0 seconds	0.001 seconds	0.001 seconds
Number of result points	513	515	515
Number of grid points	501	501	501
Number of accepted steps	361	370	348
Number of rejected steps	17	17	20
Number of f-evaluations (dynamics)	504	520	502
Number of non-linear iteration	482	494	476
Number of non-linear convergence failures	0	0	0
Number of Jacobian-evaluations	10	11	11
Number of crossing function evaluations	506	892	870
Number of model time events	6	6	6
Number of state events	0	1	1
Number of step events	0	0	0
Maximum integration order	5	5	5

The additional state events occurs due to the output onOffMea, which indicates if the device is on or off.