Olá Alessandro e demais colegas.
Você pode usar o optim para encontrar a solução aproximada. Acredito que haja uma solução mais simplificada, mas essa resolve o problema e pode ser generalizada para outras situações:
# Define os valores dos parametros que voce tem como conhecidos
b0=1
b1=1
b2=1
b3=1
b4=1
b5=1
s=1
id1=1
ab1=1
# Define expressao como funcao de x
der=function(x) {
exp(b0 + b1/s + b2/x + b3 * (id1/x) * log(ab1) + b4 * (1 - id1/x) +b5 * (1 - id1/x) * s) * (b4 * (id1/x^2) - (b3 * (id1/x^2) *log(ab1) + b2/x^2) + b5 * (id1/x^2) * s)/x - exp(b0 + b1/s +b2/x + b3 * (id1/x) * log(ab1) + b4 * (1 - id1/x) + b5 *(1 - id1/x) * s)/x^2
}
# Caso queira fazer o gráfico
curve(der(x))
# Função para otimizar. 1 é o valor inicial, e o fnscale=-1 serve para ele maximizar
optim(1, der, method="BFGS", control=list(fnscale=-1))