=====
Example 5: 2nd-order ODE - Parametric method with post-fine-tuning
=====

 

The pyCC library also supports an advanced optimization workflow referred to as post-fine-tuning. In this workflow, the code performs a two-step internal procedure. First, it optimizes the predefined model using the entire stream of reference database, exactly as shown earlier in this tutorial. Then, by simply adding extra parameters to the training routine, the code triggers an additional stage where it uses the obtained model to run forward integrations, generating simulated trajectories that are compared to the reference data. Based on this comparison, the code performs a fine-tuning of the model parameters and integrates forward again, thus performing an iterative post-fine-tuning procedure.

.. code-block:: python 


   import pycc
   import numpy as np
   import pandas as pd
   import matplotlib.pyplot as plt
   import matplotlib.ticker as ticker
   from scipy.integrate import solve_ivp
   from scipy.signal import savgol_filter
    
   # Integration of the theoretical EDO
   # We use here the method='Theoretical' using the pyCC library,
   # but we could also use standard python integrators,
   #
   # define parameters
   alpha=1.0; beta=0.2; delta=0.1; mu=0.5; omega=1.0; A=1.0
   t_span  = (0, 20)
   t_eval  = np.linspace(*t_span, 1000)
   x0     = 0.0; v0     = 0.0
   y0      = [x0, v0]  # initial conditions
   # Theoretical functions
   def F1(x_dot):
       return delta * x_dot + mu * np.tanh(500*x_dot)
   def F2(x):
       return alpha * x + beta * x**3
   def F_ext(t):
       return A * np.cos(omega * t)
   
   
   # In the following:
   # 1) we define the equations we want to integrate,
   #    (functions must be defined by f_i and params as a_i)
   # 2) we define parameters for the integrator,
   # 3) we call to pycc to make the theoretical simulation
   equations = [
       'x1_dot = x2',
       'x2_dot = F_ext - f1(x2) - f2(x1)'
   ]
   params_th = {
       # we internally use ivp_solve python function with params defined in the following 4 lines
       't_span': t_span,
       'y0': y0,
       't_eval': t_eval,
       'method': 'LSODA',
       # we need to tell to pycc how to find external/local functions
       'local_funcs': {'f1': lambda x: F1(x),'f2': lambda x: F2(x),'F_ext': lambda t: F_ext(t)},
       #'scalar_params': {'a1':2.0}
   }
   # now use pycc to integrate the EDOs
   sol,derivatives = pycc.simulate(equations,method="Theoretical", params=params_th)
   # Extract results
   time_data   = sol.t # sol is the OdeSolution returned from ivp_solve
   x1_data     = sol.y[0] # x1=x(t)
   x2_data     = sol.y[1] # x2=x_dot(t)
   F_ext_val   = F_ext(time_data)
   x1_dot_data = derivatives[0] #x1_dot=x_dot(t)
   x2_dot_data = derivatives[1] #x2_dot=x_ddot(t)
   F1_th = F1(x2_data)
   F2_th = F2(x1_data)
   
   fig, axes = plt.subplots(1, 2, figsize=(12, 6))
   axes[0].plot(time_data, x1_data, label='$x_1$ (Position)', color='blue', linewidth=2)
   axes[0].plot(time_data, x1_dot_data, label='$\\dot{x}_1=x_2$ (Velocity)', color='red', linestyle='--', linewidth=2)
   axes[0].set_title('Position and Velocity vs. Time', fontsize=16)
   axes[0].set_xlabel('Time', fontsize=12)
   axes[0].set_ylabel('Value', fontsize=12)
   axes[0].legend()
   axes[0].grid(True, which='both', linestyle='--', linewidth=0.5)
   axes[1].plot(time_data, x2_dot_data, label='$\\dot{x}_2$ (State Derivative)', color='orange', linestyle='--', linewidth=2)
   axes[1].plot(time_data, F_ext_val, label='$F_{ext}$ (External Force)', color='purple', linestyle=':', linewidth=2)
   axes[1].set_title('$x_2$, $\\dot{x}_2$ and External Force vs. Time', fontsize=16)
   axes[1].set_xlabel('Time', fontsize=12)
   axes[1].set_ylabel('Value', fontsize=12)
   axes[1].legend()
   axes[1].grid(True, which='both', linestyle='--', linewidth=0.5)
   plt.tight_layout()
   plt.show()
   
   # plotting theoretical functions f1 and f2
   fig, ax = plt.subplots(1, 2, figsize=(12, 6))
   ax[0].plot(x1_data, F1(x1_data), '.',color='orange', label="$f_1$ theory")
   ax[0].set_xlabel('$x_2$')
   ax[0].set_ylabel('$f_1(x_2)$')
   ax[0].legend()
   ax[1].plot(x2_data, F2(x2_data), '.',color='orange', label="$f_2$ theory")
   ax[1].set_xlabel('$x_1$')
   ax[1].set_ylabel('$f_2(x_1)$')
   ax[1].legend()
   plt.tight_layout()
   plt.show()
   
   
   #  Generate database that will be used for training
   # Here, we need to define all the variables that appear in eqs variable
   df = pd.DataFrame({
       #'t': time_data, # this can be added if eqs include t as an explicit variable but in this example, we encapsulated all the time dependence within F_ext_val
       'x1':x1_data,
       'x2':x2_data,
       'x1_dot':x1_dot_data,
       'x2_dot':x2_dot_data,
       'F_ext': F_ext_val
   })
   
   # Define the proposed equation, in this approach, we will use the parametric
   # approach where we want to find the ai parameters
   eqs = [
       'x1_dot = x2',
       'x2_dot = F_ext - a1 * x2 - a2* tanh(a3*x2)- a4 * x1 - a5 * x1^3',
   ]
   
   # Using Polynomial method (Poly-CC)
   # This method uses Polynomial expansions for the CCs
   # First we set the parameters
   params_poly={
     'learning_rate': 1e-3,
     'initial_param_guess':[
         {'a3':100},
     ],
     'N_order': 40,
     'n_iter': 4000,
     'eq_weights': [1.0 , 1.0],
     ## post-fine-tuning parameters are in the following:
     'fitting_forw_sim':True,  # Enables the option of refitting the coefficients using forward integr
     'n_iter_outer': 100,      # Max number of function evaluations for the outer loop
     'outer_tol': 1e-6,        # Tolerance for termination in the outer loop
     'params_simul': [
       {'local_funcs': {'F_ext': lambda t: F_ext(t)}},
       {'t_span':t_span},
       {'y0': y0},
       {'t_eval': t_eval},
       {'method': 'LSODA'},
     ],
   }
   # And here, we train the model
   models, evals , scalar_coefs = pycc.train(df=df,equations=eqs,method='Poly',params=params_poly)
   # plot obtained functions (characteristic curves CCs)
   #x_f1_cc, f1_cc, x_f2_cc, f2_cc = evals
   print(scalar_coefs)
   
   x_f1_cc = np.linspace(np.min(x1_data), np.max(x1_data), 200)
   x_f2_cc = np.linspace(np.min(x2_data), np.max(x2_data), 200)
   
   fig, ax = plt.subplots(1, 2, figsize=(12, 6))
   ax[0].plot(x_f1_cc, scalar_coefs['a1']*x_f1_cc+scalar_coefs['a2']*np.tanh(scalar_coefs['a3']*x_f1_cc), label='$f_1$ learned Poly-CC')
   ax[0].plot(x_f1_cc, F1(x_f1_cc), '--', label="$f_1$ theory")
   ax[0].set_xlabel('$x_2$')
   ax[0].set_ylabel('$f_1(x_2)$')
   ax[0].legend()
   ax[1].plot(x_f2_cc, scalar_coefs['a4']*x_f2_cc+scalar_coefs['a5']*x_f2_cc**3, label='$f_2$ learned Poly-CC')
   ax[1].plot(x_f2_cc, F2(x_f2_cc), '--', label="$f_2$ theory")
   ax[1].set_xlabel('$x_1$')
   ax[1].set_ylabel('$f_2(x_1)$')
   ax[1].legend()
   plt.tight_layout()
   plt.show()
   
   params_poly_simul = {
       'models': models, # how we enter evals instead of models
       #'evals':evals, # fi functions obtained from Poly-CC to do interpolation
       'scaling': True,
       'obtained_coefs': scalar_coefs,
       'scalar_params': scalar_coefs,
       'local_funcs': {'F_ext': lambda t: F_ext(t)},
       't_span':t_span,
       'y0': y0,
       't_eval': t_eval,
       'method': 'LSODA',
   }
   
   sol,_ = pycc.simulate(eqs, method='Poly', params=params_poly_simul)
   print(sol)
   time_sim=sol.t
   x1=sol.y[0]
   x2=sol.y[1]
   
   # Plot solution
   plt.figure()
   plt.plot(time_data, x1_data, label="x(t) th")
   plt.plot(time_sim, x1,'--', label="x(t) simulated Parametric model using Poly and post fine-tuning")
   plt.xlabel('t')
   plt.ylabel('x(t)')
   plt.legend()
   plt.show()
   

.. image:: _static/fig_ex5_1.png
   :alt: fig1
   :width: 75%
   :align: center  

.. image:: _static/fig_ex5_2.png
   :alt: fig2
   :width: 75%
   :align: center  

.. image:: _static/fig_ex5_3.png
   :alt: fig3
   :width: 75%
   :align: center  

.. image:: _static/fig_ex5_4.png
   :alt: fig4
   :width: 75%
   :align: center