I an not a maintainer of the optim package so I do not know the algorithm in detail.

But what one would expect is calculating the numerical derivative (gradient/jacobian) as first step from the starting point x0 (and other interim iterative values) to find (most) suitable direction(s) to for the (next) better solution vector iteration.

To do this, one typically calculates (f(x0 + h) - f(x0)) / h. h is very small but cannot be too small for numerical reasons, so often something like 1E-8.
So whilst the solution should follow the inequality constraints, the calculation of the derivates will lead to slightly negative values in the example provided with a very small initial starting point and therefore cause the issue.
(Note: The current linear constraints, A.' * parametervector + B >= 0 does not avoid zero in line 32/33 --> log(x/0 = Inf which is also another potential issue)


Numerically the user residual function with the log part has a discontinuity at zero (too) close to the solution space making it numerically very difficult to solve and unstable: something to avoid in optimization which prefers continuous and even better continuously differentiable functions. There will be tons of articles (tips, pitfalls, approaches) about numerical optimization of discontinuous functions in the internet.

In this specific case making the residual function more robust dealing with negative (input) values and making it at least continuous around zero/very small values will help the numerical algorithm finding a solution.
For instance a very easy solution avoiding negative p(2) is below. Now it is robust and continuous (here constant) below !e-10 ...

vsd = @(p, x) p(1)*(Tj_max+273.15)*86.1733e-6*log(x/max(1E-10, p(2))+1)+p(3)*x;

I think that the error message provided is sufficiently indicative to understand the issue (which is not an octave bug):

pars_diode = nonlin_curvefit(vsd,init,is_curve,vsd_curve,settings)
error: weighted residuals are not real
error: called from
    _lm_svd_ at line 473 column 11
    _nonlin_residmin_ at line 454 column 25
    nonlin_curvefit at line 86 column 22

The function to get the residuals is provided by the user:
vsd = @(p, x) p(1)*(Tj_max+273.15)*86.1733e-6*log(x/p(2)+1)+p(3)*x;

Numerical optimizing in general works having a starting point (init 0) and varying the (vector) points by (usually) small distances in all directions to find a better solution.

The residuals become complex according to the error message, so obviously the log in the vsd function becomes negative. Or the initial value (2) = p(2) is VERY small (9e-8 !) in the not working example. So p(2) (obviously) becomes (slightly) negative as part of the process searching for a minimum varying p.

Possible solution for this specific case is for example setting bounds / linear inequality constraints for p avoiding negative values or making the user supplied function for robust catering for negative values.




Note with regards to "it would be good to give some clarification or instruction to the user to indicate how to avoid the error by choosing a better initial guess or a different function":
a) Octave is a scientific programming language software and is for free. It assumes an user with a (reasonable) understanding of numerical analysis/programming.
There are many (scientific) books and articles about numerical analysis/optimization/programming freely available in the internet.
b) Octave as a programming language also cannot support the user with information on HOW to model user specific problems (here obviously an electronical problem) including choosing a (residual) function to optimize or its initial guess/starting value or problem specific limits/solution bounds.