first principles cooling tower improvement; (11be5caa) · Commits · Smart Buildings / EngineeringScienceBuilding

src/notebooks/Models-v3-Modelica.ipynb

0 → 100644

+484 −0

Original line number	Diff line number	Diff line
		%% Cell type:markdown id:445ce76e tags:

		# Modelica

		%% Cell type:code id:c5ea9a90 tags:

		``` python
		%reload_ext autoreload
		%autoreload 2
		import os, sys
		src_root = os.path.abspath('..')
		if src_root not in sys.path:
		sys.path.append(src_root)
		import pyfmi
		from pyfmi import load_fmu
		from systems.fmu_wrapper import FMUEnv
		import matplotlib.pyplot as plt
		import numpy as np
		```

		%% Cell type:code id:e6a2271c tags:

		``` python
		def get_variables(model: pyfmi.fmi.FMUModelCS2, match_str: str=None):
		real, integer, boolean = model.get_model_time_varying_value_references(
		filter=match_str
		)
		v = []
		v.append([])
		for r in real:
		try:
		v[-1].append((r, model.get_variable_by_valueref(r)))
		except:
		pass
		v.append([])
		for r in integer:
		try:
		v[-1].append((r, model.get_variable_by_valueref(r)))
		except:
		pass
		v.append([])
		for r in boolean:
		try:
		v[-1].append((r, model.get_variable_by_valueref(r)))
		except:
		pass
		return v

		def get_parameters(model: pyfmi.fmi.FMUModelCS2):
		# model.ge
		pass

		def get_inputs(model: pyfmi.fmi.FMUModelCS2):
		inputs = model.get_input_list()
		return tuple(inputs.keys())

		def set_parameters(model: pyfmi.fmi.FMUModelCS2, **params):
		for k, v in params.items():
		model.set(k, v)

		def get_state_vars(model: pyfmi.fmi.FMUModelCS2):
		var = get_variables(model)
		refs = []
		for category in var:
		for v in category:
		refs.append(v[1])

		refs += model.get_states_list()
		return list(set(refs))
		```

		%% Cell type:markdown id:b251a0ad tags:

		## Box

		%% Cell type:code id:81e6760c tags:

		``` python
		model = load_fmu('../systems/output/Box.fmu')
		```

		%% Cell type:code id:ed0061fa tags:

		``` python
		res = model.simulate(final_time=60, input=('F', lambda t: np.sin(t)))
		```

		%% Cell type:code id:7e7a5b09 tags:

		``` python
		res.keys()
		```

		%% Cell type:code id:79a90722 tags:

		``` python
		idx_time = res.get_column('time')
		idx_x = res.get_column('x')
		idx_vx = res.get_column('vx')
		idx_F = res.get_column('F')
		```

		%% Cell type:code id:e3e97248 tags:

		``` python
		res.data_matrix.shape
		```

		%% Cell type:code id:5321ab78 tags:

		``` python
		l, = plt.plot(res.data_matrix[idx_time], res.data_matrix[idx_F], label='F', c='r')
		plt.twinx()
		plt.plot(res.data_matrix[idx_time], res.data_matrix[idx_x], label='x', c='b')
		plt.legend(handles=[l] + plt.gca().lines)
		```

		%% Cell type:markdown id:2903c2dd tags:

		## Newton's Cooling

		%% Cell type:code id:592faf85 tags:

		``` python
		model = load_fmu('../systems/output/NewtonCooling.fmu')
		```

		%% Cell type:code id:e919c15d tags:

		``` python
		get_variables(model)
		```

		%% Cell type:code id:78ecbfba tags:

		``` python
		model.reset()
		set_parameters(model, A=1., h=1e-3, m=1.,c_p=1., T0=293.15)
		```

		%% Cell type:code id:82c32965 tags:

		``` python
		i = lambda t: 273.15 + 30 * np.sin(np.pi * t / 3600)
		```

		%% Cell type:code id:c873d279 tags:

		``` python
		res = model.simulate(start_time=0, final_time=3600, input=('T_inf', i))
		```

		%% Cell type:code id:7fa11580 tags:

		``` python
		time = res.data_matrix[res.get_column('time')]
		T_inf = res.data_matrix[res.get_column('T_inf')]
		T = res.data_matrix[res.get_column('T')]

		plt.plot(time, T_inf, ls=':', label='Ambient')
		plt.plot(time, T, ls='-', label='Output')
		plt.legend()
		```

		%% Cell type:markdown id:1d5b2986 tags:

		## PID

		%% Cell type:code id:dc9f73e5 tags:

		``` python
		model = load_fmu('../systems/output/PIDModelica.fmu', log_level=4)
		```

		%% Cell type:code id:060e4186 tags:

		``` python
		res = model.simulate(0, 10)
		plt.plot(res['time'], res['controller.y'], label='y')
		plt.plot(res['time'], res['controller.D.x'], label='D.x')
		plt.plot(res['time'], res['controller.I.y'], label='I.y')
		plt.legend()
		```

		%% Cell type:code id:e7410638 tags:

		``` python
		get_variables(model)
		```

		%% Cell type:code id:bd91ae0c tags:

		``` python
		model.set('initialOutput', 0)
		```

		%% Cell type:code id:f11be795 tags:

		``` python
		def get_pid_state(model):
		return model.get(('controller.D.x', 'controller.I.y'))
		def set_pid_state(model, derivative, integral, output=0):
		model.set('initialDerivative', derivative)
		model.set('initialIntegral', integral)
		# model.set('initialOutput', output)
		model.reset()
		print(get_pid_state(model))
		set_pid_state(model, 0, 1, 0)
		get_pid_state(model)
		```

		%% Cell type:code id:86e0588e tags:

		``` python
		model.get_log()
		```

		%% Cell type:code id:c44b6d04 tags:

		``` python
		model.reset()
		state = get_pid_state(model)
		measurement = 0
		reference = np.sin(np.linspace(0, 100, num=100) * np.pi / 50) + 1
		D = []
		I = []
		Y = []
		for i in range(100):
		model.reset()
		model.set('measurement', measurement)
		model.set('reference', reference[i])
		model.set('Ti', 1.)
		set_pid_state(model, *state)
		res = model.simulate(0, 10, options={'silent_mode': True})
		state = get_pid_state(model)
		D.append(state[0][0])
		I.append(state[1][0])
		Y.append(model.get('controller.y'))
		plt.plot(D, label='der')
		plt.plot(I, label='int')
		plt.plot(Y, label='control')
		plt.legend()
		```

		%% Cell type:markdown id:e6ea99a6 tags:

		## Cooling Tower - First principles

		%% Cell type:code id:7f9f729e tags:

		``` python
		model = load_fmu('../systems/output/Cooling.fmu')
		```

		%% Cell type:code id:327ac6d8 tags:

		``` python
		model.reset()
		model.set('ChillerHeatFlowrate', 10000)
		model.set('MassFlowrate', 5)
		res = model.simulate(start_time=0, final_time=500)
		```

		%% Cell type:code id:c4ab58d7 tags:

		``` python
		model.get_log()
		```

		%% Cell type:code id:15ab6d15 tags:

		``` python
		plt.plot(res['time'], res['TLeaving'], label='TLeaving')
		plt.plot(res['time'], res['TEntering'], label='TEntering')
		plt.plot(res['time'], res['TSetpoint'], label='TSetpoint')
		plt.plot(res['time'], res['TWetbulb'], label='TWetbulb')
		plt.legend()
		plt.twinx()
		plt.plot(res['time'], res['controller.y'], label='Controller', ls=':')
		```

		%% Cell type:code id:60febb4a tags:

		``` python
		refs = FMUEnv.get_state_vars(model)
		model.reset()
		N = 100
		step = 300
		T_wb = 290 + 5 * np.sin(np.linspace(0, N, num=N, endpoint=False) * 2 * np.pi / N)
		T_sp = T_wb + 3 * np.sin(np.linspace(0, N, num=N, endpoint=False) * 2 * np.pi / N)
		# Q_chiller =
		temps_in, temps_out = [], []
		pid_y = []
		pid_i = []
		pid_p = []
		TempCondIn, TempCondOut = 303, 303.
		# ctrl_state = get_pid_state(model)
		refs_vals = model.get(refs)
		for i in range(N):
		model.reset()
		model.set(refs, refs_vals)
		model.set('MassFlowrate', 5)
		model.set('TempCondOut', TempCondOut)
		model.set('TempCondIn', TempCondIn)
		model.set('TWetbulb', T_wb[i])
		model.set('TSetpoint', T_sp[i])
		model.set('ChillerHeatFlowrate', 50000)
		# set_pid_state(model, *ctrl_state)
		temps_in.append(TempCondOut)
		temps_out.append(TempCondIn)
		pid_y.append(model.get('controller.y')[0])
		pid_i.append(model.get('controller.I.y')[0])
		pid_p.append(model.get('controller.P.y')[0])

		res = model.simulate(start_time=0, final_time=step, options={'silent_mode': True})

		TempCondIn = res.final('TLeaving')
		TempCondOut = res.final('TEntering')
		# ctrl_state = get_pid_state(model)
		refs_vals = model.get(refs)
		```

		%% Cell type:code id:3a656a07 tags:

		``` python
		model.get_log()
		```

		%% Cell type:code id:f70bde25 tags:

		``` python
		times = np.arange(N) * step / 60
		plt.plot(times, T_sp, ls=':', label='Setpoint /K', c='g')
		plt.xlabel('Minutes')
		plt.plot(times, temps_in, label='Tower Temp In')
		plt.plot(times, temps_out, label='Tower Temp Out')
		plt.plot(times, T_wb, label='TWetbulb')
		plt.ylabel('Temperature /K)')
		plt.legend()
		plt.twinx()
		plt.plot(times, pid_y, c='r', label='PID-y')
		# plt.plot(times, pid_y, c='g', label='PID-i', ls='-.')
		# plt.plot(times, pid_y, c='b', label='PID-p', ls='-.')
		plt.legend()
		```

		%% Cell type:code id:ac6a0bfa tags:

		``` python
		from systems.fmu_wrapper import CoolingTowerFirst

		env = CoolingTowerFirst(
		model=load_fmu('../systems/output/Cooling.fmu'),
		tstep=300,
		)
		env.reset()
		env.step([290, 0.5])
		```

		%% Cell type:code id:1fd9f898 tags:

		``` python
		env.model.get('TWetbulb')
		```

		%% Cell type:code id:3d35781a tags:

		``` python
		for i in range(10):
		state = env.step([290, 15])
		print(state)
		```

		%% Cell type:markdown id:d20ec155 tags:

		## CoolingTower via Buildings lib

		`CT_FA_WB.mo`

		%% Cell type:code id:5165a20b tags:

		``` python
		model = load_fmu('../systems/output/CoolingTowerActual.fmu', log_level=4)
		```

		%% Cell type:code id:2e94f2e0 tags:

		``` python
		model.reset()
		model.set('ChillerHeatFlowrate', 10000)
		model.set('MassFlowrate', 5)
		res = model.simulate(start_time=0, final_time=500)
		```

		%% Cell type:code id:e4ab295d tags:

		``` python
		plt.plot(res['time'], res['tow.TLvg'])
		plt.plot(res['time'], res['tow.TAir'])
		plt.twinx()
		plt.plot(res['time'], res['tow.y'], ls=':', c='r')
		```

		%% Cell type:code id:ccfe02c0 tags:

		``` python
		model.get_log()
		```

		%% Cell type:code id:e73941b6 tags:

		``` python
		refs = FMUEnv.get_state_vars(model)
		model.reset()
		N = 100
		step = 300
		T_wb = 290 + 5 * np.sin(np.linspace(0, N, num=N, endpoint=False) * 2 * np.pi / N)
		T_sp = T_wb + 3 * np.sin(np.linspace(0, N, num=N, endpoint=False) * 2 * np.pi / N)
		temps_in, temps_out = [], []
		pid_y = []
		pid_i = []
		pid_p = []
		T0 = 303.
		# ctrl_state = get_pid_state(model)
		refs_vals = model.get(refs)
		for i in range(N):
		model.reset()
		model.set(refs, refs_vals)

		model.set('TempCondOut', T0)
		model.set('MassFlowrate', 0.5)
		model.set('TSetpoint', T_sp[i])
		model.set('TWetbulb', T_wb[i])
		model.set('ChillerHeatFlowrate', 50000)
		# set_pid_state(model, *ctrl_state)

		res = model.simulate(start_time=0, final_time=step, options={'silent_mode': True})
		refs_vals = model.get(refs)
		# ctrl_state = get_pid_state(model)
		temps_in.append(T0)
		temps_out.append(res.final('tow.TLvg'))
		T0 = model.get('vol.T')[0] # tower entering temperature
		pid_y.append(model.get('controller.y')[0])
		pid_i.append(model.get('controller.I.y')[0])
		pid_p.append(model.get('controller.P.y')[0])
		```

		%% Cell type:code id:5919fb3a tags:

		``` python
		times = np.arange(N) * step / 60
		plt.xlabel('Minutes')
		l, =plt.plot(times, T_sp, ls=':', label='TSetpoint', c='g')
		plt.plot(times, temps_in, label='Tower Temp In')
		plt.plot(times, temps_out, label='Tower Temp Out')
		plt.plot(times, T_wb, label='T_wetbulb')
		plt.ylabel('Temperature /K)')
		lines = plt.gca().lines
		plt.twinx()
		plt.plot(times, pid_y, c='r', label='PID-y')
		plt.plot(times, pid_y, c='g', label='PID-i', ls='-.')
		plt.plot(times, pid_y, c='b', label='PID-p', ls='-.')
		plt.legend(handles=plt.gca().lines + lines)
		```

		%% Cell type:code id:25c9128f tags:

		``` python
		from systems.fmu_wrapper import CoolingTowerActual

		env = CoolingTowerActual(
		model=load_fmu('../systems/output/CoolingTowerActual.fmu'),
		tstep=300,
		)
		env.reset()
		```

		%% Cell type:code id:a2033e23 tags:

		``` python
		for i in range(10):
		state = env.step([290, 0.5])
		print(state)
		```

src/notebooks/Models-v3.ipynb

deleted100644 → 0

+0 −24

Original line number	Diff line number	Diff line
		%% Cell type:code id:appropriate-ribbon tags:

		``` python
		import notebook_setup

		import numpy as np
		import matplotlib.pyplot as plt
		import pandas as pd
		```

		%% Cell type:code id:contained-assumption tags:

		``` python
		class X:
		def __init__(self):
		pass
		def __getitem__(self, key):
		print(key)
		X()[1]
		```

		%% Output

		1

src/systems/CT_FA_WB.mo

+9 −11

Original line number	Diff line number	Diff line
		@@ -9,8 +9,9 @@ partial model PartialStaticTwoPortCoolingTower

		parameter Modelica.SIunits.MassFlowRate m_flow_nominal = 0.5
		"Design water flow rate";
		parameter Modelica.SIunits.HeatFlowRate q_flow_evaporator = 0.5m_flow_nominal4200*5;
		parameter Modelica.SIunits.HeatFlowRate q_flow_evaporator = m_flow_nominal*4184;
		parameter Real volume = 0.5;
		parameter Real T_start=300;

		replaceable Buildings.Fluid.HeatExchangers.CoolingTowers.BaseClasses.CoolingTower
		tow
		@@ -52,6 +53,7 @@ partial model PartialStaticTwoPortCoolingTower
		redeclare package Medium = Medium_W,
		m_flow_nominal=m_flow_nominal,
		V=volume, // volume /m^3 of mixing volume/water loop
		T_start=T_start, // Starting temperature of water
		energyDynamics=Modelica.Fluid.Types.Dynamics.FixedInitial);

		Buildings.Fluid.Sources.Boundary_pT exp(
		@@ -91,11 +93,12 @@ model CoolingTowerActual

		parameter Real TSetpoint = 273.15 + 18;
		parameter Real TWetbulb = 273.15 + 15;
		parameter Real T0 = 300;
		parameter Real TempCondIn = 300;
		parameter Real TempCondOut = 300;

		parameter Real MassFlowrate = 0.5;
		parameter Real Volume = 0.5;
		parameter Real ChillerHeatFlowrate = 0;
		parameter Real Volume = 0.5;

		parameter Real k = 1;
		parameter Real Ti = 10;
		@@ -108,24 +111,19 @@ model CoolingTowerActual
		m_flow_nominal=MassFlowrate
		),
		q_flow_evaporator=ChillerHeatFlowrate,
		volume=Volume
		volume=Volume,
		T_start=TempCondOut
		);

		// Modelica.Blocks.Interfaces.RealInput TAir(
		// final min=0,
		// final unit="K",
		// displayUnit="degC")
		// "Entering air wet bulb temperature";

		Modelica.Blocks.Continuous.LimPID controller(
		k=k,
		Ti=Ti,
		Td=Td,
		// reverseActing=false,
		initType=Modelica.Blocks.Types.InitPID.InitialState,
		xd_start=initialDerivative, xi_start=initialIntegral,
		yMin=0, yMax=1)
		"Controller for tower fan";

		initial equation

		equation

src/systems/CoolingTower.mo

+38 −13

Original line number	Diff line number	Diff line
		@@ -18,25 +18,40 @@ end NewtonCooling;
		model ThermalCapacitance "A model of thermal capacitance"
		parameter Modelica.SIunits.HeatCapacity C "Thermal capacitance = specific capacity x mass";
		parameter Real HeatFlowRate=0 "Heat flowing into the system";
		Modelica.Thermal.HeatTransfer.Interfaces.HeatPort body;
		Modelica.Thermal.HeatTransfer.Interfaces.HeatPort tower;
		equation
		C*der(body.T) = body.Q_flow + HeatFlowRate;
		C*der(body.T) = tower.Q_flow;
		end ThermalCapacitance;


		model EvaporativeCooling "Cooling to ambient conditions given airspeed and temperature"
		model Tower "Cooling to ambient conditions given airspeed and temperature"
		import Connectors = Modelica.Blocks.Interfaces;
		parameter Real coolingConstant=500 "Constant of proportionality J/(K s)";
		parameter Real capacitance;
		parameter Real T_wb=300 "Wetbulb temperature";

		Connectors.RealInput airspeed "Speed of air over water";
		Modelica.Thermal.HeatTransfer.Interfaces.HeatPort load;
		Modelica.Thermal.HeatTransfer.Interfaces.HeatPort load_in, load_out;

		equation
		// TODO: Equation governing evaporation under airflow and water flow
		load.Q_flow = coolingConstant * (1 + airspeed) * (load.T - T_wb);

		end EvaporativeCooling;
		load_in.Q_flow + load_out.Q_flow = coolingConstant * (1 + airspeed) * (load_in.T - T_wb);
		capacitance * der(load_out.T) = load_in.Q_flow + load_out.Q_flow;
		end Tower;


		model Chiller "Exchange with chiller"
		import Connectors = Modelica.Blocks.Interfaces;

		parameter Real HeatFlowrate=0;
		Modelica.Thermal.HeatTransfer.Interfaces.HeatPort load_in, load_out;
		parameter Modelica.SIunits.HeatCapacity capacitance "Thermal capacitance = specific capacity x mass";

		equation
		load_out.Q_flow + load_in.Q_flow = -HeatFlowrate;
		capacitance * der(load_out.T) = HeatFlowrate;
		end Chiller;


		model Fan "Cooling Tower Fan"
		@@ -57,7 +72,8 @@ class Cooling

		parameter Real TSetpoint = 273.15 + 18;
		parameter Real TWetbulb = 273.15 + 15;
		parameter Real T0=300 "initial water temperature entering tower";
		parameter Real TempCondIn=300 "initial water temperature exiting tower";
		parameter Real TempCondOut=300 "initial water temperature entering tower";

		parameter Real MassFlowrate = 0.5;
		parameter Real Volume = 0.5;
		@@ -69,9 +85,10 @@ class Cooling
		parameter Real Td = 0;
		parameter Real initialIntegral=0, initialDerivative=0;

		EvaporativeCooling evap(coolingConstant=coolingConstant, T_wb=TWetbulb);
		Tower tow(coolingConstant=coolingConstant, T_wb=TWetbulb, capacitance=Volume * 1e3 * 4200);
		Chiller chi(capacitance=Volume * 1e3 * 4200, HeatFlowrate=ChillerHeatFlowrate);
		// Capacitance = mass * specific heat capacity = volume * density * c_p
		ThermalCapacitance water(C=Volume * 1e3 * 4200, HeatFlowRate=ChillerHeatFlowrate);
		// ThermalCapacitance water(C=Volume * 1e3 * 4200, HeatFlowRate=ChillerHeatFlowrate);

		LimPID controller(
		k=k,
		@@ -83,10 +100,18 @@ class Cooling
		) "Controller for tower fan";

		initial equation
		water.body.T = T0;
		// water.body.T = T0;
		tow.load_in.T = TempCondOut;
		tow.load_out.T = TempCondIn;
		// chi.load_in.T = TempCondIn;
		// chi.load_out.T = TempCondOut;
		equation
		evap.airspeed = controller.y;
		controller.u_m = -water.body.T; // Negative sign to make PID reverse acting
		controller.u_m = -tow.load_out.T; // Negative sign to make PID reverse acting
		controller.u_s = -TSetpoint;
		connect(water.body, evap.load);
		tow.airspeed = controller.y;
		// chi.HeatFlowrate = ChillerHeatFlowrate;
		// connect(water.body, tow.load);
		connect(chi.load_out, tow.load_in);
		connect(tow.load_out, chi.load_in);
		// connect(tow.load_out, tow.load_in);
		end Cooling;

src/systems/CoolingTower_single.mo

0 → 100644

+74 −0

File added.

Preview size limit exceeded, changes collapsed.

Admin message