Electrochemical-process/main.py at main · AndrosovAS/Electrochemical-process · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
import numpy as np

import intvalpy as ip
ip.precision.extendedPrecisionQ = False

from data import TMES, WMES


################################################################################
################################################################################
# Введём модельную функцию у которой будем восстанавливать параметры
# Вычислим её градиенты
# W(.., t) = (x0 + x1*t) / (x2 + x3*t)

# +----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+
# model = 0
model = lambda a, x: (x[0] + x[1]*a) / (x[2] + a)

dx0 = lambda a, x: 1 / (x[2] + a)
dx1 = lambda a, x: a / (x[2] + a)
dx2 = lambda a, x: -(x[0] + x[1]*a) / (x[2] + a)**2

# +----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+
# model = 1
# model = lambda a, x: (x[0] + K1*x[1]*a) / (x[2] + x[1]*a)
#
# dx0 = lambda a, x: 1 / (x[2] + x[1]*a)
# dx1 = lambda a, x: a* (K1*x[2] - x[0]) / (x[2] + x[1]*a)**2
# dx2 = lambda a, x: -(x[0] + K1*x[1]*a) / (x[2] + x[1]*a)**2

# +----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+----+

grad = np.array([dx0, dx1, dx2])
# количество неизвестных
N = len(grad)


################################################################################
# Поскольку все входные измерения различны (кроме финального) и независимы,
# то можно объединить в одну большую выборку

tmes = np.array([mes for exp in TMES for mes in exp])
wmes = np.array([mes for exp in WMES for mes in exp])

# Известно, что данные неточные
eps_t = 0
eps_w = 0.0255

tmes = tmes + ip.Interval(-eps_t, eps_t)
wmes = wmes + ip.Interval(-eps_w, eps_w)

# Известна связь между коэффициентами
K1 = 0.75


################################################################################
# Алгоритмически вычислим выбросы
# Для этого воспользуемся методом распознающего функционала
# Поскольку погрешность по времени отсутствует, то можно взять Tol

# Сначала нам нужно найти точку, которая в наибольшей степени согласуется с данными
# Интересно сравнить с МНК-точкой

# Для нахождения такой точки проведём глобальную оптимизацию и найдём глобальный
# максимум распознающего функционала
# В случае дробно-линейной функции функционал Tol является квазивогнутой функцией,
# поэтому риска зациклиться в локальном экстремуме нет

# Следующим шагом смотрим на образующие функционала
# Если их значение сильно отклоняется от большинства, то можно утверждать,
# что это является выбросом
# Это можно определить с помощью инквартильного метода или с помощью метода стандартного отклонения
# Возможно, некоторые выбросы скрываются за более явными выбросами, поэтому
# может возникнуть необходимость несколько раз повторить алгоритм.

# Введём функции для количественного выявление выбросов
def interquartile(data):
    q25, q75 = np.percentile(data, 25), np.percentile(data, 75)
    iqr = q75 - q25
    cut_off = iqr * 1.5

    lower, upper = q25 - cut_off, q75 + cut_off
    return np.argwhere((data < lower) | (data > upper)).flatten()

def standard_deviations(data):
    # Set upper and lower limit to 3 standard deviation
    std, mean = np.std(data), np.mean(data)
    cut_off = std * 3

    lower, upper = mean - cut_off, mean + cut_off
    return np.argwhere((data < lower) | (data > upper)).flatten()

# Задаём начальное приближение
x0 = np.array([171.32052537,   1.00037412, 224.51843273])

# Повторяем описанный алгоритм пока не иссякнет максимальное количество итераций
# или массив выбросов не станет пустым
nit = 0
maxiter = 10

while nit < maxiter:
    success, x, f = ip.nonlinear.Tol(model, tmes, wmes, maxQ=True, grad=grad, x0=x0, tol=1e-8, maxiter=5000)
    print(success)

    tt = wmes.rad - (wmes.mid - model(tmes, x)).mag
    outliers_index = standard_deviations(tt)

    if len(outliers_index) > 0:
        index = np.delete(np.arange(tmes.shape[0]), outliers_index)
        tmes, wmes = tmes[index], wmes[index]
        print('outliers_index: ', outliers_index)
    else:
        break
    nit += 1

# Таким образом получена чистая выборка
# Посмотрим на вектор x и максимальное значение функционала
print('x: ', x)
print('f: ', f)