Ajuste de mínimos quadrados: scipy.optimize.curve_fit

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. :)

Deixe um comentário