Files
MSE-PI-E2EEDA-Plein-de-eeee…/model/Study_CO2_concentration_window_opening.py
Alison Lecointre a84ad89e02 created models
2026-06-04 12:38:11 +02:00

293 lines
11 KiB
Python

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 26 16:19:28 2026
@author: alisonlecointre
"""
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.pyplot as ticker
import numpy as np
import glob
import os
# %%=============================
# File loading
# ===============================
def load_file(file_name):
sheet= pd.read_excel(file_name, sheet_name=None) #loading file, sheat_name=None : return dictionary containing df for each sheet
for sheet_name, df in sheet.items(): #loading all sheets from the file
df = df.dropna(axis=1, how="all") #remove empty column "all" : all elements missing
df.columns = df.columns.str.strip() #remove leading character
sheet[sheet_name]=df
return sheet
def load_all_files(folder): #load all .xlsx files from the folder
files = glob.glob(f"{folder}/*.xlsx")
all_data={} #create dictionnary of all dictionnary of each file
for file in files:
name = os.path.basename(file)
all_data[name] = load_file(file)
return all_data
folder = os.getcwd()
data = load_all_files(folder)
# %%=============================
# Parameters Provence classrooms
# ===============================
df_Provence_rooms = pd.DataFrame({
"classrooms": ["A2","A3","A4","A5","A6","A7"],
"volume": [308, 480, 326,274,323,272],
"n_student_max": [40, 78, 48,32,58,36],
"windows_surface": [3.23, 22.12, 3.23, 22.12, 3.23, 5.53]
})
print (df_Provence_rooms)
# %%=============================
# Calculation
# ===============================
# ===============================
# Air volume per person
# ===============================
def calculation_air_volum_per_pers (df) :
air_volume_per_pers = df["Mean room volume (m3)"]/ df["Mean number of students (-)"]
return air_volume_per_pers
calculation = {} #new dictionnary to store the air volum per person
for file_name, sheet in data.items() : #for each file including sheets from this file and all included in data
calculation [file_name] = {}
for sheet_name, df in sheet.items() : #df in all sheets
calculation[file_name][sheet_name] = calculation_air_volum_per_pers(df)
# ================================
# New data frame for simulation model of CO2 concentration increase over time
# ================================
CO2_prod_th = 0.0155/60 #m3/min <=> ajusted at 15.5 l/h after crossing data with simaria model [20 l/h CO2 production per person during normal activity (cf. article 16 822.113 Ordonnance 3 du 18 août 1993 relative à la loi sur le travail (OLT 3) (Protection de la santé))
def model_no_windows_opening (CO2_t0=None, N=None, V=None, df=None) :
t = np.arange(0, 181, 1)
if df is not None :
CO2_t0 = df["CO2 rate (ppm)"].iloc[0]
N = df["Mean number of students (-)"].iloc[0]
V = df["Mean room volume (m3)"].iloc[0]
CO2_t_model = CO2_t0 + ((N * CO2_prod_th * t *1e6) / V)
return t, CO2_t_model
model = {}
for file_name, sheet in data.items() : #for each file including sheets from this file and all included in data
model [file_name] = {}
for sheet_name, df in sheet.items() : #df in all sheets
t, CO2_t_model = model_no_windows_opening(df=df)
model [file_name][sheet_name] = {"t": t, "CO2_t_model": CO2_t_model}
# ================================
# Find threshold time
# ================================
threshold = 1400
def find_threshold_time(t, CO2_t_model, threshold):
indices = np.where(CO2_t_model >= threshold)[0]
if len(indices) == 0:
return None # never reached
return t[indices[0]]
threshold_time = {}
for file_name, sheet in model.items():
threshold_time[file_name] = {}
for sheet_name, values in sheet.items():
t = values["t"]
CO2_t_model = values["CO2_t_model"]
threshold_time[file_name][sheet_name] = find_threshold_time(t, CO2_t_model, threshold)
# ================================
# New data frame for simulation model of CO2 concentration decrease over time (using SIA 382/1 norm)
# ================================
v_air = 0.5 #m/s Assumption (cf. article)
def model_windows_opening_CO2 (S_windows=None, volume=None, df=None):
t = np.arange(0, 1, 0.01)
if df is not None:
for i, row in df.iterrows():
S_windows = df_Provence_rooms["windows_surface"].iloc[i]
volume = df_Provence_rooms["volume"].iloc[i]
q_air_change = S_windows * v_air * 3600 #m3/h
CO2_ext = 400 #ppm
CO2_int_0 = threshold #ppm
CO2_t_model_ventilation = ((CO2_int_0 - CO2_ext - ((0.001 * (CO2_prod_th * 1000*60))/q_air_change))* np.exp((-q_air_change * t)/volume) + CO2_ext + ((0.001 * (CO2_prod_th * 1000*60))/q_air_change) ) #ppm CO2_t_model_ventilation = ((CO2_int_0 - CO2_ext - ((0.001 * (CO2_prod_thx))/q_air_change))* np.exp((-q_air_change * t)/volume) + CO2_ext + ((0.001 * (CO2_prod_thx))/q_air_change) ) #ppm
t = t*60
return t, CO2_t_model_ventilation
model_with_windows = {}
for i, row in df.iterrows():
t, CO2_t_model_ventilation = model_windows_opening_CO2(df=row.to_frame().T)
model_with_windows[i] = pd.DataFrame({"CO2_t_model_ventilation" : CO2_t_model_ventilation, "t" : t})
# ================================
# Find window opening duration
# ================================
def find_windows_opening_time (t,CO2_t_model_ventilation):
windows_opening_duration = np.interp(405, CO2_t_model_ventilation[::-1] , t[::-1])
return windows_opening_duration
# ================================
# Air quality thresholds for graphs
# ================================
thresholds = {
"Good air quality":(0, 1400,"green"),
"Bad air quality": (1400, 2000,"orange"),
"Really bad air quality": (2000, np.inf,"red"),
}
def plot_air_quality_zones (ax, thresholds):
xmin, xmax = ax.get_xlim()
ymin, ymax = ax.get_ylim()
for label, (y_min, y_max, color) in thresholds.items():
if y_max != np.inf:
ax.axhline(y=y_max, linestyle="--", color=color)
y_top = min(y_max, ymax)
y_bottom = max(y_min, ymin)
y_mid = (y_bottom + y_top) / 2
ax.axhspan(y_bottom, y_top, color=color, alpha=0.2)
ax.text(
(xmin + xmax) / 2,
y_mid,
label,
ha="center",
va="center",
color="black",
bbox=dict(facecolor="white", alpha=0.5)
)
# %% =============================
# Visualisation
# ================================
# ================================
# Graph time evolution of CO2 concentration open data
# ================================
for file_name, sheet in data.items():
fig0, ax = plt.subplots(figsize=(10, 7))
for sheet_name, df in sheet.items():
t = model[file_name][sheet_name]["t"]
CO2_t_model = model[file_name][sheet_name]["CO2_t_model"]
t_threshold = threshold_time[file_name][sheet_name]
ax.plot(df["time (min)"], df["CO2 rate (ppm)"], label=f"{file_name} - {sheet_name}")
ax.plot(t,CO2_t_model, label = f"{file_name} - {sheet_name} simulated", linestyle="--")
if t_threshold is not None :
ax.axvline(x = t_threshold, label = "time reaching threshold", linestyle=":", color="red")
ax.legend(
loc="upper center",
bbox_to_anchor=(0.45, -0.08),
ncol=3 if len(ax.get_legend_handles_labels()[1]) > 1 else 1
)
fig0.tight_layout()
ax.set_ylim(0, 5000)
plot_air_quality_zones (ax, thresholds)
ax.grid(True, alpha=0.7)
ax.set_xlabel("time (min)")
ax.set_ylabel("CO2 concentration (ppm)")
ax.set_title("CO2 concentration over time")
# ================================
# Graph time evolution of CO2 concentration after window opening
# ================================
fig1, ax= plt.subplots(figsize=(10, 6))
for i in model_with_windows:
t = model_with_windows [i]["t"]
CO2_t_model_ventilation = model_with_windows[i]["CO2_t_model_ventilation"]
ax.plot(t, CO2_t_model_ventilation, label=df_Provence_rooms.loc[i, "classrooms"])
ax.legend()
ax.xaxis.set_major_locator(ticker.MultipleLocator(2))
ax.set_xlim(0, 20)
ax.grid(True, alpha=0.7)
ax.set_xlabel("time (min)")
ax.set_ylabel("CO2 concentration (ppm)")
ax.set_title("CO2 concentration over time after window opening - with open data")
# %%==========================================================================================
# User
# ==========================================================================================
print("\n")
name_classroom = input("Classroom number (A2 to A7) : ")
volume_room = df_Provence_rooms.loc[df_Provence_rooms["classrooms"] == name_classroom, "volume"].values
S_windows_room = df_Provence_rooms.loc[df_Provence_rooms["classrooms"] == name_classroom, "windows_surface"].values
n_student = float(input("Number of students : "))
initial_CO2_level = float(input("Initial CO2 level (ppm): "))
while initial_CO2_level < 400:
print("The value must be >= 400 ppm, which correspond to the normal CO2 concentration in air. Try again.")
initial_CO2_level = float(input("Initial CO2 level (ppm): "))
#Create CO2 concentration evolution over time using simulation model and give time reaching threshold
t_user, CO2_t_model_user = model_no_windows_opening(CO2_t0=initial_CO2_level, N=n_student, V=volume_room, df=None)
t_threshold = find_threshold_time(t_user, CO2_t_model_user, threshold)
t_user2, CO2_t_model_ventilation_user = model_windows_opening_CO2(S_windows =S_windows_room, volume = volume_room, df=None)
t_user3 = t_user2 + t_threshold #start at the time reaching threshold
windows_opening_duration = find_windows_opening_time (t_user2, CO2_t_model_ventilation_user)
print(f"Time reaching threshold of {threshold} ppm =", round(t_threshold),"minutes")
print(f"Open windows for {windows_opening_duration:.0f} min after reaching threshold")
#Display CO2 concentration evolution over time
reponse = input("Display CO2 evolution over time? (yes/no) : ").strip().lower()
while reponse not in ["yes", "no"]:
print("Answer with yes or no")
reponse = input("Display CO2 evolution over time? (yes/no) : ").strip().lower()
if reponse == "yes":
fig10, ax = plt.subplots(figsize=(10, 6))
ax.plot(t_user, CO2_t_model_user , label="CO2 concentration over time (simulated)")
ax.plot(t_user3, CO2_t_model_ventilation_user, label="with opened windows", color="brown")
if t_threshold is not None :
ax.axvline(x = t_threshold, label = "time reaching threshold", linestyle=":", color="red")
ax.legend(
loc="upper center",
bbox_to_anchor=(0.45, -0.07),
ncol=3 if len(ax.get_legend_handles_labels()[1]) > 1 else 1
)
fig10.tight_layout()
plot_air_quality_zones(ax, thresholds)
ax.set_title("Time evolution of CO2 concentration - Simulated")
ax.set_xlabel("Time (min)")
ax.set_ylabel("CO2 concentration (ppm)")
ax.set_xlim(0, 150)
ax.set_ylim(400, 5000)
ax.grid(True)