Solving a system of a changing number of nonlinear equations

I'm trying to solve the following system of nonlinear equations for a spacecraft thermal modeling project:

This system grows in size with the number of nodes, and because of the sums, I really don't want to have to code the equations in manually (especially since I'll be altering the number of nodes as I add detail to the model--if I had a fixed set I might just suck it up and do it).

My advisor has recommended I use fsolve to work this out, but while I think he'd be right for using it for a fixed number of nodes, for a dynamic version I'm not so sure.

This is the cleanest version of what I've tried so far.

Currently, it gets through the loop, but then 'crashes' on trying to solve infinitely, and I think is getting stuck in a NaN/Infinity loop between trustnleqn.m and sym.m.

fun = @test1
x0 = T_0
X = fsolve(fun,x0)

function One = test1(X)
load 'AppxSteadyStateTestInput.mat'
boltzmann = 5.67E-8;
N = numel(nodes);
One = zeros(N, 1);

for i = 1:1:N
    for j = 1:1:N
    One(i) = -(Q_external(i)+Q_d(i) - boltzmann*eps_outside(i)*A_r(i)*X(i) - sum(Conductances(i,j)*(X(i)-X(j)))-boltzmann*sum(F_times_A(i)*epsilon_inside(i,j)*(X(i)^4-X(j)^4)));
    end

end

end

I'm not sure if I'm trying to do something impossible, or if I just haven't found the right way to do it yet.

I've also attempted the same math with vpasolve, but I get wildly incorrect numbers, but at least this code runs, so maybe that's just an input issue or some other minor problem....

load 'AppxSteadyStateTestInput.mat'
boltzmann = 5.67E-8;

% Conductances = zeros(size(Conductances))
N = numel(nodes);

x = sym('T', [N 1] );

T_0 = [273; 273; 273; 273; 273; 273; 273; 273; 293; 273]; %initial conditions

results = sym('results',[N 1]);
SUM1 = sym(zeros(N,N));
SUM2 = sym(zeros(N,N));

for i = 1:1:N
    for j = 1:1:N
        syms k
        SUM1(i,j) = symsum(Conductances(i,j)*(x(i)-x(j)), k, i, j);
        SUM2(i,j) = symsum(F_times_A(i)*epsilon_inside(i,j)*(x(i)^4-x(j)^4),k,i,j);
        results(i) = -(Q_external(i)+Q_d(i) - boltzmann*eps_outside(i)*A_r(i)*x(i) - SUM1(i,j) -boltzmann*SUM2(i,j));
    end
end

disp('coderun')

vpasolve(results,x)

[Walter Roberson]

        SUM1(i,j) = symsum(Conductances(i,j)*(x(i)-x(j)), k, i, j);

Note that symsum() is completely unable to have the variable of summation be used for indices. So your symsum() here is equvailent to

        EXPRESSION = Conductances(i,j)*(x(i)-x(j));
        SUM1(i,j) = symsum(EXPRESSION, k, i, j);

which in turn is equivalent to

       EXPRESSION = Conductances(i,j)*(x(i)-x(j));
       SUM1(i,j) = EXPRESSION .* (j-i+1);

Using symsum() as a replacement for straight multiplication is inefficient and confusing.

The same remarks apply to the second symsum -- as written it is equivalent to multiplication.

[Matt J]

Your code is difficult to interpret, since you use variables that are entirely different from your mathematical expression,

Also, it appears that the code contains an extra term,

boltzmann*eps_outside(i)*A_r(i)*X(i)

not present in the original equation.

Answer 1

[Torsten]

The solver doesn't converge. You should check parameters and equations.

load 'AppxSteadyStateTestInput.mat'
fun = @(x)test1(x,boltzmann,nodes,Q_external,Q_d,eps_outside,A_r,Conductances,F_times_A,epsilon_inside);
x0 = T_0;
X = fsolve(fun,x0,optimset('MaxFunEvals',1000000,'MaxIter',1000000))

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+06.

X = 10×1

1.0e+04 * 2.0289 2.0291 2.0291 2.0292 2.0291 2.0293 2.0291 2.0291 2.0291 2.0291

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function One = test1(X,boltzmann,nodes,Q_external,Q_d,eps_outside,A_r,Conductances,F_times_A,epsilon_inside)
  N = numel(nodes);
  One = zeros(N, 1);
  for i = 1:N
    One(i) = Q_external(i) + Q_d(i) - boltzmann*eps_outside(i)*A_r(i)*X(i);
    for j = 1:N
       One(i) = One(i) - Conductances(i,j)*(X(i)-X(j)) - boltzmann*F_times_A(i)*epsilon_inside(i,j)*(X(i)^4-X(j)^4);
    end
  end
end

Answer 2

[Matt J]

Because your equations are 4-th order polynomials, they are apt to have multiple roots. You probably won't get a physically sensible result unless you impose some bounds on T (and maybe other constraints). Here is how that might look:

load 'AppxSteadyStateTestInput.mat';

W0=-boltzmann.*eps_outside.*A_r;
W1=-Conductances;
W2=-boltzmann.*F_times_A.*epsilon_inside;
Q = Q_external+Q_d;

W1=diag(W0)+diag(sum(W1,2))-W1;
W2=diag(sum(W2,2))-W2;

fun=@(T) deal( Q + W1*T + W2*T.^4   ,   W1+4*W2.*(T.^3')  );

opts=optimoptions('lsqnonlin',MaxFunEval=Inf,MaxIterations=1e4, ...
                   OptimalityTol=1e-12, FunctionTol=1e-12,StepTol=1e-12,...
                   SpecifyObjectiveGradient=true);

lb=zeros(size(T_0));

[T, resnorm,residuals]=lsqnonlin(fun,T_0,lb+200,lb+350, opts)

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

T = 10×1

348.9869 350.0000 349.7822 349.7822 350.0000 348.9869 328.5763 336.2230 346.2899 328.5763

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

resnorm = 6.4422e+03

residuals = 10×1

25.3814 25.3814 25.3814 25.3814 25.3814 25.3814 25.3814 25.3814 25.3814 25.3814

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Nicholas]

I've spent the afternoon playing with this code working to understand it and the matrix formulation--very cool stuff! Have not heard of/used Deal before.

Unfortunately, it seems that the solver gets a residual norm in the 10^7 range, which as I understand from the lsqnonlin documentation page, is rather bad. I would think that the next logical step, then, is to apply additional constraints, but given that this set of equations is only a function of one variable vector, I can't immediately think of any additional constraints to apply--or at least, any that would enter the equation at the solver level (there are a number of constraints in the input file, but they shouldn't be changing as a function of temperature, except maybe conductance, and since I am attempting to replicate a model from literature, I'm using the paper's values exactly as inputs, which means that the conductances are fixed).

Thank you for the succinct and clear answer, nonetheless.

[Matt J]

Unfortunately, it seems that the solver gets a residual norm in the 10^7 range

That's not what I'm getting. I've edited my answer to give the resnorm output, and it is 6.4422e+03.

which as I understand from the lsqnonlin documentation page, is rather bad.

It might be bad, but I don't think anything in the doc page could have told you that. There is nothing in the magnitude of the resnorms alone that tells you the quality of the solution. You can change the resnorm, without changing the quality of the solution, just by multiplying all your equations by a constant.

[Matt J]

and since I am attempting to replicate a model from literature, I'm using the paper's values exactly as inputs

If this has already been done in previous literature, it would be wise, as a sanity-check, to plug in their solution for T and check if it solves the equations -- the equations as implemented by you.

I don't see how the equations, as you've presented them to us, have a well-behaved solution.

[Torsten]

@Nicholas

and since I am attempting to replicate a model from literature, I'm using the paper's values exactly as inputs,

Maybe you also have the solution for T from the paper, and you could evaluate your function with it to see where you might have made a mistake in the implementation.

I think that somewhere in your equations, there should appear a fixed environmental temperature in the chain of resistances to close the system.

Although your system is stationary and includes radiative heat transfer, this answer might help:

https://uk.mathworks.com/matlabcentral/answers/2176334-how-to-solve-a-thermal-model-need-help?s_tid=srchtitle