Ajustar dados experimentais usando algum modelo com função analítica ou numérica usando ajustes por mínimos quadrados.
Para ver uma aplicação experimental deste código, ver Diodo Schottky com resistências parasitas.
O código discutido aqui usa partes do código explicado no Método 2 de Solução numérica para uma equação transcendental.
Este texto não pretende explicar o método de ajustes por mínimos quadrados. Porém, se você tem interesse em saber um pouco mais sobre o método, leia Ajuste por mínimos quadrados: versão papel e caneta
Usando scipy.optimize.curve_fit
A documentação da função curve_fit: scipy.optimize.curve_fit
Para obter a versão completa do código explicado abaixo acesse o repositório desse código no GitHub. (ou baixe pacote zip).
import matplotlib.pyplot as plt import numpy as np import scipy.optimize as optm import __main__ from scipy import special from matplotlib import font_manager
No #BLOCO 1 (linha 3) inicializamos o módulo scipy.optimize que contém a função curve_fit.
def JVdiode(vv, Rs, Rp, m, i0): kT = __main__.kT log_ii_teo = [] for v in vv: r = (v-i0*Rp)/(Rp+Rs) a0 = (Rp+Rs)/(i0*Rp)*np.exp(-v/(m*kT)) c = Rs/(m*kT) z = c/a0*np.exp(-c*r) i = r+(1/c)*special.lambertw(z) if(abs(i.real) < abs(i0)): log_ii_teo.append(np.log10(abs(i.real))) else: log_ii_teo.append(np.log10(i0)) nlog_ii_teo=np.array(log_ii_teo) return nlog_ii_teo
Para usar a função curve_fit precisamos definir a função do modelo de ajuste que possa ser chamada com a variável independente como primeiro argumento e os parâmetros de ajuste em sequência. Na linha 1 do #BLOCO 2 definimos JVdiode em função do vetor vv que contem os valores de voltagem. Rs, Rp, m e i0 são parâmetros do modelo que podem ser alterados para ajustar os dados experimentais. Detalhes importantes:
- No código usado anteriormente para resolver o problema da equação transcendental (veja aqui), utilizamos matrizes 1D tipo lists do python para guardar os valores de voltagem e corrente. Para utilizar a função curve_fit faz-se necessário utilizar matrizes do tipo numpy arrays. Então todas as matrizes são criadas como listas usando a append (linhas 11 e 13) e depois transformadas em numpy arrays usando np.array() (linha 14).
- Na função a ser chamada pela curve_fit todos os argumentos, exceto o primeiro (vv neste caso), serão usados como parâmetros de ajustes e não podem ser constantes. Outras parâmetros constantes devem ser importados da rotina principal usando __main__ (linha 2).
- Neste exemplo, utilizo a função W de Lambert para resolver uma equação transcendental. A função lambertw pode, em certas situações, retornar valores complexos. Para evitar problemas de convergência na função curve_fit eu decidi usar apenas a componente real da corrente calculada (linha 11).
- Para os meus dados experimentais, os valores de corrente variam ordens de grandeza. Para que o ajuste seja feito de forma adequada em todos o range de voltagem, inclusive onde a corrente é muito baixa comparada com os valores mais altos de corrente, é preciso realizar o ajuste do logarítimo da corrente. Desta forma, o valor retornado pela função JVdiode (linha 15) é calculado usando a função log10 (linhas 11 e 13).
- Por fim, encontrei dificuldades calculando log10 de valores muito baixos de corrente. Por isso decidi limitar o uso da função ao valor de i0 (ver condicional na linha 10)
# PARAMETROS FIXOS T = 330 #in K k = 8.617e-5 #in eV/K kT = k*T # PARAMETROS VARIÁVEIS Rs = 4.0e2 # in Ohms Rp = 8.5e7 # in Ohms m = 1.35 i0 = 1.0e-10 # in Amperes # INICIALIZANDO VETORES vv = [] ii_exp = [] ii_exp_abs = [] ii_teo = [] log_ii_exp = [] # directory = "diretorio/para/seu/dado/experimental/" # Directory of files. End with / path = directory+"nome_do_arquivo.csv" #arquivo em duas colunas separadas por virgula. look = 0 with open(path, 'r') as data: for line in data.readlines(): if(look): vv.append(float(line.split(",")[0])) ii_exp.append((float(line.split(",")[1]))) ii_exp_abs.append(abs((float(line.split(",")[1])))) log_ii_exp.append(np.log10(abs(float(line.split(",")[1])))) if(line[0:7] == "Voltage"): look = 1 nvv = np.array(vv) nii_exp = np.array(ii_exp) nii_exp_abs = np.array(ii_exp_abs) nlog_ii_exp=np.array(log_ii_exp) # ESCOLHER INTERVALO DE INTERESSE EM TORNO DO CENTRO DO GRAFICO ran = 0 nvv = nvv[ran:-1-ran] nii_exp = nii_exp[ran:-1-ran] nii_exp_abs = nii_exp_abs[ran:-1-ran] nlog_ii_exp = nlog_ii_exp[ran:-1-ran]
No #BLOCO 3 não tem muito o que falar além das linhas 32 a 35 onde matrizes lists são transformadas em matrizes numpy arrays.
x0 = np.array([Rs, Rp, m, i0]) sol, cov = optm.curve_fit(JVdiode, nvv, nlog_ii_exp, x0) Rs = sol[0] Rp = sol[1] m = sol[2] i0 = sol[3] nlog_ii_teo = JVdiode(nvv, Rs, Rp, m, i0) for i in nlog_ii_teo: ii_teo.append(10.0**i)
No #BLOCO 4 é onde toda ação acontence. Na linha 1 x0 é uma numpy array com todas os parâmetros de ajustes da função JVdiode. No meu caso, inicializei os valores no #BLOCO 3. Na linha 2 chamamos a função curve_fit que por sua vez chama a função JVdiode, usa a matriz de variável independente (vv), do logarítmo da corrente como matriz de dados experimentais (nlog_ii_exp) e x0. A função curve_fit retorna duas matrizes sol e cov. A matrix sol contém os valores otimizados dos parâmetros de ajuste. A matriz cov contém informações relacionados ao erro do ajuste. Depois que o ajuste é feito, podemos calcular os valores da função ajustada recuperando os parâmetros (linhas 4 a 7), chamando a função JVdiode (linha 9) e calculando a corrente resultante (linha 12).
Encontrou algum erro? Entra em contato e ajuda a corrigir. :)