<jupyter_start><jupyter_text>Designing an optimal Coffee extraction with Python Necessary libs:<jupyter_code>import scipy.integrate as integrate
from scipy.optimize import differential_evolution
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D<jupyter_output><empty_output><jupyter_text>Solving the ODE Model data:<jupyter_code>phi_h = 0.8272
k_sv1 = 27.35e-6
k_sv2 = 322.49e-6
l_l = 282e-6
D_h = D_v = 2.2e-9
m = 30e-6
c_sat = 212.4
c_s = 1400.0
kappa = 3.1
phi_c0 = 0.143435
alpha = 0.1833
beta = 0.0447
r_s = 16.94<jupyter_output><empty_output><jupyter_text>Initial conditions:<jupyter_code>c_h0 = 0.
c_v0 = 183.43
phi_v0 = 0.6444
psi_s0 = 1<jupyter_output><empty_output><jupyter_text>EDO's system right-hand-side (RHS):<jupyter_code>def coffee(t, X, alpha, beta):
    #Variables
    c_h, c_v, phi_v, psi_s = X
    
    # Coffee's desorption parameters
    desorption_1 = alpha * phi_v ** (4. / 3.) * D_v * 6. / (k_sv2 * l_l)
    desorption_2 = beta * 12. * D_h * phi_c0 / (k_sv1 * m)
    
    # Time derivatives
    # c_h
    c_h_prime = - (1. - phi_h) * desorption_1 * (c_h - c_v) + (1 - phi_h) * desorption_2 * (c_sat - c_h) * psi_s
    c_h_prime = (1. / phi_h) * c_h_prime
    # c_v
    c_v_prime = - desorption_1 * (c_h - c_v) - (c_v / c_s) * desorption_2 * (c_sat - c_h) * psi_s
    c_v_prime = (1. / phi_v) * c_v_prime
    # psi_s
    psi_s_prime = - (r_s / c_s) * desorption_2 * (c_sat - c_h) * psi_s
    # phi_v
    phi_v_prime = - (1. / r_s) * psi_s_prime
    return c_h_prime, c_v_prime, phi_v_prime, psi_s_prime
def coffee(t, X, alpha, beta):
    #Variables
    c_h, c_v, phi_v, psi_s = X
    
    # Coffee's desorption parameters
    desorption_1 = alpha * phi_v ** (4. / 3.) * D_v * 6. / (k_sv2 * l_l)
    desorption_2 = beta * 12. * D_h * phi_c0 / (k_sv1 * m)
    
    # Time derivatives
    # c_h
    c_h_prime = - (1. - phi_h) * desorption_1 * (c_h - c_v) + (1 - phi_h) * desorption_2 * (c_sat - c_h) * psi_s
    c_h_prime = (1. / phi_h) * c_h_prime
    # psi_s
    psi_s_prime = - (r_s / c_s) * desorption_2 * (c_sat - c_h) * psi_s
    # c_v
    c_v_prime = - desorption_1 * (c_h - c_v) + (c_v / r_s) * psi_s_prime
    c_v_prime = (1. / phi_v) * c_v_prime
    # phi_v
    phi_v_prime = - (1. / r_s) * psi_s_prime
    return c_h_prime, c_v_prime, phi_v_prime, psi_s_prime<jupyter_output><empty_output><jupyter_text>Time step parameters:<jupyter_code>t0 = 0  # tempo inicial
tmax = 600. # tempo máximo<jupyter_output><empty_output><jupyter_text>Solving with RADAU IIA method:<jupyter_code>solution_RADAU = integrate.solve_ivp(
    lambda t, X: coffee(t, X, alpha, beta), 
    t_span=(t0, tmax), 
    y0=(c_h0, c_v0, phi_v0, psi_s0,), 
    method='Radau'
)

solution_RADAU<jupyter_output><empty_output><jupyter_text>Gathering the solutions:<jupyter_code>times, variables = solution_RADAU.t, solution_RADAU.y<jupyter_output><empty_output><jupyter_text>Separating the solutions:<jupyter_code>c_h, c_v, phi_v, psi_s = variables<jupyter_output><empty_output><jupyter_text>Plotting the solutions:<jupyter_code>plt.plot(times, c_h)
plt.show()<jupyter_output><empty_output><jupyter_text>Calibration: Reading the data from experimental measures:<jupyter_code>exp_data = pd.read_csv('../data/coffee_conc.csv', decimal=",")

exp_data
plt.plot(exp_data.t, exp_data.c_h, 'o-')
plt.xlabel('Time [s]')
plt.ylabel(r'$c_h$ [kg/m$^3$]')
plt.grid(True, linestyle='--')
plt.xlim([0, 500])
plt.ylim([0, 40])
plt.show()<jupyter_output><empty_output><jupyter_text>Calibration objective function:<jupyter_code>def coffee_solution(x, t_eval=[]):
    if len(x) > 2:
        raise ValueError('x must have length equal to 2: alpha and beta parameters')
        
    # Getting fitting parameters
    alpha, beta = x
    
    # Solving the ODE system
    if len(t_eval) == 0:
        raise ValueError("Please provide times to evaluate the function")
    else:
        try:
            solution_RADAU = integrate.solve_ivp(
                lambda t, X: coffee(t, X, alpha, beta), 
                t_span=(0, t_eval.max()), 
                y0=(c_h0, c_v0, phi_v0, psi_s0,), 
                method='Radau',
                t_eval=t_eval.values
            )
        except:
            return np.array([np.inf] * len(t_eval))
    
    c_h, c_v, phi_v, psi_s = solution_RADAU.y
    
    if len(c_h) != len(t_eval):
        c_h = np.array([np.inf] * len(t_eval))
    
    return c_h
def calibration_least_squares(pars, exp_data):
    residual = exp_data.c_h.values - coffee_solution(pars, t_eval=exp_data.t)
    least_squares_evaluation = np.sum(residual ** 2.0)
    return least_squares_evaluation<jupyter_output><empty_output><jupyter_text>Fitting by least square minimization:<jupyter_code>bounds = [(0, 1), (0, 1)]
result = differential_evolution(
    calibration_least_squares,
    bounds=bounds,
    args=(exp_data,),
    disp=True,
    popsize=70,
    strategy='best2exp',
    tol=1e-5,
    recombination=0.9,
    mutation=0.5,
    maxiter=50
)
alpha, beta = result.x

alpha, beta
solution_RADAU = integrate.solve_ivp(
            lambda t, X: coffee(t, X, alpha, beta), 
            t_span=(t0, tmax), 
            y0=(c_h0, c_v0, phi_v0, psi_s0,), 
            method='Radau'
        )

time = solution_RADAU.t
c_h, c_v, phi_v, psi_s = solution_RADAU.y
plt.plot(time, c_h, label='Simulated')
plt.plot(exp_data.t, exp_data.c_h, 'o-', label='Experimental')
plt.xlabel('Time [s]')
plt.ylabel(r'$c_h$ [kg/m$^3$]')
plt.grid(True, linestyle='--')
plt.xlim([0, 500])
plt.ylim([0, 40])
plt.legend()
plt.savefig('calibration_model1.png')
plt.show()
bounds = [(-alpha, alpha), (-beta, beta)]
result = differential_evolution(
    calibration_least_square,
    bounds=bounds,
    args=(exp_data,),
    disp=True,
    popsize=70,
    strategy='best2bin',
    tol=1e-5,
    recombination=0.9,
    mutation=0.5,
    maxiter=50
)
alpha, beta = result.x
solution_RADAU = integrate.solve_ivp(
            lambda t, X: coffee(t, X, alpha, beta), 
            t_span=(t0, tmax), 
            y0=(c_h0, c_v0, phi_v0, psi_s0,), 
            method='Radau'
        )

time = solution_RADAU.t
c_h, c_v, phi_v, psi_s = solution_RADAU.y
plt.plot(time, c_h, label='Simulated')
plt.plot(exp_data.t, exp_data.c_h, 'o-', label='Experimental')
plt.xlabel('Time [s]')
plt.ylabel(r'$c_h$ [kg/m$^3$]')
plt.grid(True, linestyle='--')
plt.xlim([0, 500])
plt.ylim([0, 40])
plt.legend()
plt.savefig('calibration_model1.png')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14202004626060138, 'fraction_numerical': 0.031457208943716265, 'mean_word_length': 3.385395537525355, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '179647', 'license': 'MIT', 'n_tokens_mistral': 2982, 'n_tokens_neox': 2786, 'n_words': 519, 'path': 'notebook/explanation_model1.ipynb', 'repo_name': 'volpatto/coffeebreak', 'size': '6485'}