Projekt 1

"""
    Obliczanie podstawowych wskaźników jakości regulacji dla podanego przebiegu

    @param _tp Okres próbkowania
    @param _sp Lista wartości zadanych
    @param _t Lista wartości czasu
    @param _e Lista wartości odchyłki

    @return Słownik wskaźników jakości regulacji {czas regulacji, odchyłka statyczna, ISE, IAE}
"""
def calculate_quality_indicators(_tp, _t, _sp, _e):
    # Czas regulacji (przyjmujemy tolerancję 2% wartości zadanej)
    settling_time = None
    tolerance = 0.02 * _sp  # 2% wartości zadanej na końcu symulacji
    stable_region_start = None
    # Szukamy momentu, od którego sygnał pozostaje w tolerancji
    for i in range(len(_t) - 1, -1, -1):
        if abs(_e[i]) > tolerance:
            stable_region_start = i + 1  # Czas regulacji jest po ostatnim przekroczeniu tolerancji
            break
    if stable_region_start is not None and stable_region_start < len(_t):
        settling_time = _t[stable_region_start]
    else:
        settling_time = math.inf
    # Odchyłka statyczna (różnica między wartością zadaną a wyjściową w stanie ustalonym)
    # Przyjmujemy ostatnią wartość jako stan ustalony
    static_error = _e[-1]
    # Całka z kwadratu odchyłki (ISE - Integral Squared Error)
    ise = sum([e ** 2 * _tp for e in _e])
    # Całka z wartości bezwzględnej odchyłki (IAE - Integral Absolute Error)
    iae = sum([abs(e) * _tp for e in _e])
    return {
        "st": settling_time,
        "se": static_error,
        "ISE": ise,
        "IAE": iae
    }
    

"""
    Wyświetla przebiegi i wskaźniki podsumowujące jakośc regulacji:
        - wykres SP i PV
        - wykres e
        - wykres CV

    @param _t Lista wartości czasu
    @param _sp Lista wartości zadanych
    @param _pv Lista wartości wielkości regulowanej
    @param _e Lista wartości odchyłki
    @param _cv Lista wartości sygnału sterującego
    @param _eval Wynik evaluacji z funkcji calculate_quality_indicators(). Domyślnie: None
    @param _elim Dodatkowe ograniczenie progowe do narysowania na przebiegu e. Domyślnie: None
    @param _splim Lista granic do narysowania na przebiegu SP. Domyślnie: None
    @param _pv_bis Lista drugich wartości wielkości PV, np. odtwaranych z modelu
"""
def show_evaluation(_t, _sp, _pv, _e, _cv, _eval=None, _elim=None, _splim=None, _pv_bis=None):
    plt.figure()
    fig, ax = plt.subplots(nrows=3, ncols=1, figsize=(9, 15))

    SP_H = [x + 0.02 for x in _sp]
    SP_L = [x - 0.02 for x in _sp]

    limSP_H = [x * 0.02 for x in _sp]
    limSP_L = [x * (-0.02) for x in _sp]

    ax[0].plot(_t, _sp, label='SP', color='r')
    ax[0].plot(_t, _pv, label='PV', color='g')
    if _pv_bis is not None:
        ax[0].plot(_t, _pv_bis, label='PV_pred', color='palegreen')
    ax[0].plot(_t, SP_H, linestyle='', marker='.', markersize=2, label='2%', color='b', )
    ax[0].plot(_t, SP_L, linestyle='', marker='.', markersize=2, color='b')
    ax[0].legend()
    if _eval is not None:
        ax[0].set_title(f"st: {_eval['st']:.2f}")
    if _splim is not None:
        for pos, v in enumerate(_splim):
            ax[0].plot([0, _t[-1]], [_splim[pos], _splim[pos]], label=('lim' if pos==0 else ''), color='y')

    ax[1].plot(_t, _e, label='e', color='g')
    ax[1].plot(_t, limSP_H, linestyle='', marker='.', markersize=2, label='2%', color='b')
    ax[1].plot(_t, limSP_L, linestyle='', marker='.', markersize=2, color='b')
    if _elim is not None:
        ax[1].plot([0, _t[-1]], [_elim, _elim], label='lim', color='y')
        ax[1].plot([0, _t[-1]], [-_elim, -_elim], color='y')
    ax[1].legend()
    if _eval is not None:
        ax[1].set_title(f"ISE: {_eval['ISE']:.3f}")

    ax[2].plot(_t, _cv, label='CV', color='g')
    ax[2].legend()
    

def fitness_func(ga_instance, solution, solution_idx):
    # Symulacja
    system.in_var_val('kp', solution[0])
    proc_vars = system.simulate(tend, 'control')
    # fe = filter(dt, proc_vars['e'])
    # Wyznaczenie wskaźników oceny jakości regulacji
    eval = calculate_quality_indicators(dt, proc_vars['t'], sp, proc_vars['e'])
    fitness1 = 1.0 / (eval['st'] + 0.000001)
    fitness2 = 1.0 / (eval['ISE'] + 0.000001)
    progress_bar(solution_idx+1, sol_per_pop)
    print(f" Step: {solution_idx+1} Solution: kp={solution[0]} Eval - st: {eval['st']}, ISE: {eval['ISE']}")
    return [fitness1, fitness2]
    

proc_vars = system.simulate(tend, 'control')
evaluation = calculate_quality_indicators(dt, proc_vars['t'], sp, proc_vars['e'])
show_evaluation(proc_vars['t'], proc_vars['SP'], proc_vars['PV'], proc_vars['e'], proc_vars['CV'], evaluation, tend)
    

# Parametry AG
num_generations = 100
num_parents_mating = 10
sol_per_pop = 20
#   dostrajać będziemy tylko jeden parametr: kp
num_genes = 1
init_range_low = 0.1
init_range_high = 100

ga_instance = pygad.GA(num_generations=num_generations,
                       num_parents_mating=num_parents_mating,
                       sol_per_pop=sol_per_pop,
                       num_genes=num_genes,
                       mutation_type = None,
                       fitness_func=fitness_func,
                       parent_selection_type='nsga2',
                       init_range_low=init_range_low,
                       init_range_high=init_range_high)
ga_instance.run()