Calculation of parameters

[tlmerbecks]

Ciao Francesco, thanks for merging my pull request for #587!

I have been thinking about this issue a bit further, and have gotten the feeling that such errors could in future be avoided by devolving the calculation of parameters from the residual and post-processing calculation.

For example, in this case, the error arose because the calculation of ttd_min had been replicated in various places and a typo had crept into one of these instances. Instead, we could define a function that calculates the parameter, which can then be called whenever needed, i.e. for the residual and post-processing calculations.

This could be implemented by adding another attribute to the parameter data container, namely a function to calculate the parameter, which is then called by the residual function. In the example below, "resi" (currently "func") is the function to calculate the residual and "func" is the function to calculate the parameter.

def get_parameters():
    return {'eta_s': dc_cp(
                            min_val=0, max_val=1, num_eq=1,
                            deriv=self.eta_s_deriv,
                            resi=self.eta_s_resi,
                            func=self.eta_s_func, 
                            latex=self.eta_s_func_doc)
                }

This could potentially also streamline the post-processing calculations, by allowing all properties to be calculated by looping over the parameters and calling the respective "func" function.

Hope this makes sense

[fwitte]

Good morning @tlmerbecks,

happy new year and thank you for the follow-up.

I have also thought about something into this direction (but not as nice as you suggested). I really like it and I think I will implement it together with a couple of other changes:

• Rename deriv to derivative (maybe even jacobian, what do you think?) and func to residual (as proposed by you). Writing down the complete names sounds nicer somehow.
• Add the a new field calculate_value, calculate or function (what do you like best?) that will calculate the actual value.
• Implement a method in the Component class that will return the specified parameter value minus the value of the function of the previous bullet. This will be the default residual function, e.g. default_residual(), but not be utilized for all methods.

It will look something like this:

def default_residual(self, parameter, **kwargs):
    return parameter.val - parameter.function(**kwargs)

And the parameter settings like this:

@staticmethod
def get_parameters():
    return {
        'eta_s': dc_cp(
            min_val=0, max_val=1, num_eq=1,
            derivative=self.eta_s_deriv,
            residual=self.eta_s_resi,
            function=self.calculate_eta_s, 
        ),
        'P': dc_cp(
            min_val=0, max_val=1, num_eq=1,
            derivative=self.energy_balance_deriv,
            residual=self.default_residual,
            function=self.calculate_P, 
        )
    }

The reason why one parameter would NOT use the default function can be of three reasons:

First, the function might not actually be formulated in the way value - calculate_value() for improved performance, i.e. eta_s_func of a compressor. The partial derivative to h_2 in this formulation requires more resources.

$$0 = \eta_s - \frac{h_{2,s}-h_1}{h_2-h_1}$$

In this formulation it is very simple, i.e. just $\eta_s$:

$$0 = \eta_s \dot (h_2-h_1) - (h_{2,s}-h_1)$$

tespy/src/tespy/components/turbomachinery/compressor.py

Lines 192 to 219 in 8d47499

	def eta_s_func(self):
	r"""
	Equation for given isentropic efficiency of a compressor.
	Returns
	-------
	residual : float
	Residual value of equation.
	.. math::
	0 = -\left( h_{out} - h_{in} \right) \cdot \eta_{s} +
	\left( h_{out,s} - h_{in} \right)
	"""
	i = self.inl[0]
	o = self.outl[0]
	return (
	(o.h.val_SI - i.h.val_SI) * self.eta_s.val - (
	isentropic(
	i.p.val_SI,
	i.h.val_SI,
	o.p.val_SI,
	i.fluid_data,
	i.mixing_rule,
	T0=None
	) - self.inl[0].h.val_SI
	)
	)

tespy/src/tespy/components/turbomachinery/compressor.py

Lines 240 to 262 in 8d47499

	def eta_s_deriv(self, increment_filter, k):
	r"""
	Partial derivatives for isentropic efficiency.
	Parameters
	----------
	increment_filter : ndarray
	Matrix for filtering non-changing variables.
	k : int
	Position of derivatives in Jacobian matrix (k-th equation).
	"""
	i = self.inl[0]
	o = self.outl[0]
	f = self.eta_s_func
	if self.is_variable(i.p, increment_filter):
	self.jacobian[k, i.p.J_col] = self.numeric_deriv(f, 'p', i)
	if self.is_variable(o.p, increment_filter):
	self.jacobian[k, o.p.J_col] = self.numeric_deriv(f, 'p', o)
	if self.is_variable(i.h, increment_filter):
	self.jacobian[k, i.h.J_col] = self.numeric_deriv(f, 'h', i)
	if self.is_variable(o.h, increment_filter):
	self.jacobian[k, o.h.J_col] = self.eta_s.val

The second reason is when convergence stability (avoiding mathematical errors) is required, i.e. in kA functions of heat exchanger:

tespy/src/tespy/components/heat_exchangers/base.py

Lines 443 to 495 in 8d47499

	def calculate_td_log(self):
	i1 = self.inl[0]
	i2 = self.inl[1]
	o1 = self.outl[0]
	o2 = self.outl[1]
	# temperature value manipulation for convergence stability
	T_i1 = i1.calc_T()
	T_i2 = i2.calc_T()
	T_o1 = o1.calc_T()
	T_o2 = o2.calc_T()
	if T_i1 <= T_o2:
	T_i1 = T_o2 + 0.01
	if T_i1 <= T_o2:
	T_o2 = T_i1 - 0.01
	if T_i1 <= T_o2:
	T_o1 = T_i2 + 0.02
	if T_o1 <= T_i2:
	T_i2 = T_o1 - 0.02
	ttd_u = T_i1 - T_o2
	ttd_l = T_o1 - T_i2
	if ttd_u == ttd_l:
	td_log = ttd_l
	else:
	td_log = (ttd_l - ttd_u) / math.log((ttd_l) / (ttd_u))
	return td_log
	def kA_func(self):
	r"""
	Calculate heat transfer from heat transfer coefficient.
	Returns
	-------
	residual : float
	Residual value of equation.
	.. math::
	0 = \dot{m}_{in,1} \cdot \left( h_{out,1} - h_{in,1}\right) +
	kA \cdot \frac{T_{out,1} -
	T_{in,2} - T_{in,1} + T_{out,2}}
	{\ln{\frac{T_{out,1} - T_{in,2}}{T_{in,1} - T_{out,2}}}}
	"""
	return (
	self.inl[0].m.val_SI * (
	self.outl[0].h.val_SI - self.inl[0].h.val_SI
	) + self.kA.val * self.calculate_td_log()
	)

The calculate_td_log function, checks if there are impossible values. During iterations there might be impossible values (negative ttd_l or ttd_u) which would cause the log to raise an error. In order to avoid that, the values are manipulated in the process. At the end of the calculation in the postprocessing, a different method is used to actually calculate the td_log and kA which will prompt a warning, if any of the parameters remained physically infeasible until the end of the calculation.

The third reason would be a parameter, which does not utilize a value like we are discussion here, but is just an expression (e.g. the Stodola law for turbine):

tespy/src/tespy/components/turbomachinery/turbine.py

Lines 249 to 277 in 8d47499

	def cone_func(self):
	r"""
	Equation for stodolas cone law.
	Returns
	-------
	residual : float
	Residual value of equation.
	.. math::
	0 = \frac{\dot{m}_{in,ref} \cdot p_{in}}{p_{in,ref}} \cdot
	\sqrt{\frac{p_{in,ref} \cdot v_{in}}{p_{in} \cdot v_{in,ref}}}
	\cdot \sqrt{\frac{1 - \left(\frac{p_{out}}{p_{in}} \right)^{2}}
	{1 - \left(\frac{p_{out,ref}}{p_{in,ref}} \right)^{2}}} -
	\dot{m}_{in}
	"""
	n = 1
	i = self.inl[0]
	o = self.outl[0]
	vol = i.calc_vol(T0=i.T.val_SI)
	return (
	- i.m.val_SI + i.m.design * i.p.val_SI / i.p.design
	* (i.p.design * i.vol.design / (i.p.val_SI * vol)) ** 0.5
	* abs(
	(1 - (o.p.val_SI / i.p.val_SI) ** ((n + 1) / n))
	/ (1 - (self.pr.design) ** ((n + 1) / n))
	) ** 0.5
	)

So the envisioned solution would be:

• Every parameter of components (and we can transfer the same logic to connections and possibly other things for the future) uses:
  • the derivative function as is
  • the default_residual function in case it does not fall below any of the examples above
  • the calculate_value function, in case a value can be calculated
• The solving then calls the residual functions of used equations, which (in case of the default_residual) function will be using the calculate_value functions.
• In postprocessing, the calculate_value functions of every parameter that has been assigned one ... if calculate_value is not None ... will be executed and written into the parameter value itself.

Maybe you can give me a quick feedback of what you think about the solution :).

Thank you!

[tlmerbecks]

Ciao @fwitte,

Happy New Year to you too!

My only concern with the approach outlined above is that the renaming of the parameters (e.g. func to residual or deriv to jacobian) could create backwards-compatibility issues...?

Personally, I perfer renaming func to residual and deriv to jacobian, as it aligns the naming of these quantities to the naming already used at the network level.

As for the new calculate function, I am undecided... part of me is thinking that the field should just be called value instead of calculate_value as all of the other fields also calculate something; however, value is somewhat non-descript and also mirrors val and val_SI... Besides the others you suggested what do you think of calc_value or definition.

Other than that, I fully agree with the solution outlined in the bullet points in you previous post.

[fwitte]

My only concern with the approach outlined above is that the renaming of the parameters (e.g. func to residual or deriv to jacobian) could create backwards-compatibility issues...?

It certainly would, but these will come at some point no matter what you do, thus the change will be a new major version, 0.8.x. this will also be the case with component property units, I want to implement that soon! Same for the concept with the variable space reduction...

Personally, I perfer renaming func to residual and deriv to jacobian, as it aligns the naming of these quantities to the naming already used at the network level.

I will stick with them then.

As for the new calculate function, I am undecided... part of me is thinking that the field should just be called value instead of calculate_value as all of the other fields also calculate something; however, value is somewhat non-descript and also mirrors val and val_SI... Besides the others you suggested what do you think of calc_value or definition.

Since func and deriv do not use verbs, I would then go with definition I think.

Thank you, have a nice day