how to implement reduced order model in MATLAB

[tyler davis]

Hello - I have been working with Evgenii to reduce a sample thermal
finite element model. I have the reduced system matrices

The thermal model is linear (for simplicity) and first order, of the
form:

E dT/dt + A T = B

I have approximately 20,000 nodes (degrees of freedom)

the reduced model is of the form:

E_r dz/dt + A_r z = B_r
T_out = C_r z

which now has only 30 degrees of freedom. (The _r subscript denotes
the reduced matrix size)

I have 6 output nodes of interest. Two of them represent the location
of Dirichlet (defined-temperature) boundary conditions. One of the
BCs is fixed with time, and one of them changes throughout the
simulation.

I have a simple MATLAB script that performs the time integration of
the reduced matrices using the ode15s solver. The initial state vector
is set to the initial temperature of the thermal nodes.

Following simulation, the resulting state vector z is multiplied by
the transformation matrix C_r to yield T_out for the 6 nodes of
interest.

I can't see how to implement the two boundary conditions during the
simulation of the reduced model. The Dirichlet fixed temperature on
one of the output nodes is 100°C, but the initial temperatures are all
25°C. I don't understand how to insert this BC into the time-
integration simulation of the reduced model.

Hopefully this explains my problem and thanks for any help you can
offer

[Luis Miguel Silveira]

Hi,

I really do not know how you generated the reduced model. I presume it was with some sort of projection. If this is the case, then there is a projector, V (a tall thin matrix, few long columns) that relates T with z_r, such as

T = V * z

If the procedure was such that V has orthogonal columns (as is common in many MOR algorithms), then V^T V = I which implies that

z_r = V^T * T

Thus, if you have initial conditions for T, you can translate that into z_r using V^T.

Not sure if this helps. Best regards,

-L. Miguel Silveira

[Evgenii Rudnyi]

> I can't see how to implement the two boundary conditions during the
> simulation of the reduced model. The Dirichlet fixed temperature on
> one of the output nodes is 100�C, but the initial temperatures are all
> 25�C. I don't understand how to insert this BC into the time-
> integration simulation of the reduced model.

First I would recommend you to make ambient of 0 degree. This
corresponds to a temperature transformation

T_new = T - 25

and then you add 25 after simulation of the reduced model. This is the
easiest way to deal with not zero initial conditions. Otherwise you
should do what Miguel has suggested but this is in the general case
actually tricky. In this respect you may want to look at

Model Order Reduction for Scanning Electrochemical Microscope: The
treatment of nonzero initial conditions
http://modelreduction.com/doc/papers/feng04IEEEsensors.pdf

As for the not zero Dirichlet, this should be considered as an extra
input. You first make Dirichlet zero and apply power dissipation. This
will give you the first input. Then you make power dissipation zero and
apply not zero Dirichlet equal to one. This gives you the second input
that you can scale afterwards in the reduced model.

It does not make sense to define Dirichlets as outputs, as they are known.

[tyler davis]

thank you for the post

The model was reduced using MOR 2.51

I believe I am implementing a similar algorithm to compute the initial
state vector that corresponds to the initial Temperature vector.
However using this method, the initial Temperature vector is still 0
degrees, rather than 25 degrees.

The dimension of the transformation vector C_r is [6 x 30]. This
corresponds to 6 output nodes of interest and a reduced state=space
model with 30 degrees of freedom.

Perhaps it would help to post the MATLAB code I am using. It is very
simple.

____________________________________________________________________________________
%this file exercises the reduced model generated by 'MOR for ANSYS'

%clear old and define global matrices
clear A B C E
global A B E num_nodes

%read in matrices
[A]=mmread('mor.A');
[B]=mmread('mor.B');
[C]=mmread('mor.C');
[E]=mmread('mor.E');

%determine matrix dimensions
num_nodes = size(B,1);
num_outputs =size(C,1);
tinitial=0;tfinal=10800; %set integration range

%set intial conditions
T=zeros(num_outputs,1); %initialize temperature vector
T(:)=25; %set intial temperature
z0=T/C; %calculate initial state vector corresponding to T_initial
vector

%run ODE solution
atol=ones(num_nodes,1)*1e-5;
options = odeset('RelTol',1e-3,'AbsTol',atol);
[t,z] = ode15s(@firstorderODE,[tinitial tfinal],z0);

%plot figure
T=C*transpose(z); %turn down the output dimension by ode15 call
figure(1);
plot(t,T)
xlabel('time (s)');
ylabel('temperature (°C)')
title('mor30 transient results');
____________________________________________________________________________________________
function dzdt = firstorderODE(t,z)
global A B E num_nodes
dzdt = zeros(num_nodes,1);
dzdt = (B-(A*z))\E;
dzdt=transpose(dzdt);
______________________________________________________________________________________________

[tyler davis]

>
> As for the not zero Dirichlet, this should be considered as an extra
> input. You first make Dirichlet zero and apply power dissipation. This
> will give you the first input. Then you make power dissipation zero and
> apply not zero Dirichlet equal to one. This gives you the second input
> that you can scale afterwards in the reduced model.
>
> It does not make sense to define Dirichlets as outputs, as they are known.

I agree that Dirichlet nodes are not outputs to the reduced model, but
rather should be inputs. The temperature of the Dirichlet nodes
affects all the other temperatures in the reduced model, since they
are linked through conductors. I don't understand your method, but I
don't think it will work for my model. I want to use the reduced
model to evaluate different settings of the Dirichlet nodes (that
represent different environments). I literally want to "plug in" the
defined temperatures of the boundary nodes, and have them propagate
through the other non-boundary nodes in the reduced model.

Perhaps it would help to show you a sample of the results from the
full-fidelity model (click on the link for image)
http://picasaweb.google.com/lh/photo/N4c8gDdQAj7NsisvbyjbxQ?feat=directlink

As you can see, the other nodes in the model respond to changes in the
temperature of the Dirichlet boundaries. I cannot use a simple before/
after offset to simulate this effect. It has to be in the model at
runtime.

On Mar 31, 1:41 pm, Evgenii Rudnyi <use...@rudnyi.ru> wrote:
> > I can't see how to implement the two boundary conditions during the
> > simulation of the reduced model.  The Dirichlet fixed temperature on
> > one of the output nodes is 100 C, but the initial temperatures are all
> > 25 C.  I don't understand how to insert this BC into the time-
> > integration simulation of the reduced model.
>
> First I would recommend you to make ambient of 0 degree. This
> corresponds to a temperature transformation
>
> T_new = T - 25
>
> and then you add 25 after simulation of the reduced model. This is the
> easiest way to deal with not zero initial conditions. Otherwise you
> should do what Miguel has suggested but this is in the general case
> actually tricky. In this respect you may want to look at
>
> Model Order Reduction for Scanning Electrochemical Microscope: The

> treatment of nonzero initial conditionshttp://modelreduction.com/doc/papers/feng04IEEEsensors.pdf

[Evgenii Rudnyi]

>> As for the not zero Dirichlet, this should be considered as an extra
>> input. You first make Dirichlet zero and apply power dissipation. This
>> will give you the first input. Then you make power dissipation zero and
>> apply not zero Dirichlet equal to one. This gives you the second input
>> that you can scale afterwards in the reduced model.
>>
>> It does not make sense to define Dirichlets as outputs, as they are known.
>
> I agree that Dirichlet nodes are not outputs to the reduced model, but
> rather should be inputs. The temperature of the Dirichlet nodes
> affects all the other temperatures in the reduced model, since they
> are linked through conductors. I don't understand your method, but I
> don't think it will work for my model. I want to use the reduced
> model to evaluate different settings of the Dirichlet nodes (that
> represent different environments). I literally want to "plug in" the
> defined temperatures of the boundary nodes, and have them propagate
> through the other non-boundary nodes in the reduced model.
>
> Perhaps it would help to show you a sample of the results from the
> full-fidelity model (click on the link for image)
> http://picasaweb.google.com/lh/photo/N4c8gDdQAj7NsisvbyjbxQ?feat=directlink
>
> As you can see, the other nodes in the model respond to changes in the
> temperature of the Dirichlet boundaries. I cannot use a simple before/
> after offset to simulate this effect. It has to be in the model at
> runtime.

Let us write some equations. First we sort degrees of freedom, xn -
normal, xd - Dirichlets. Then we have

(Kn Knd) (xn) = F
(Kdn Kd) (xd) = F_reactive

When you fix Dirichlets, then you have reactive forces (heat flow) on
the Dirichlet DoFs but we are not interested in them now. What is
important though is to understand that in the second equation xd are
known and F_reactive are unknown. This is a big difference with the the
first equation. Now you need just remove xd from the first equation
because they are known. So we have

Kn xn = F - Knd*xd

As I have said, the not zero Dirichlet acts as the extra load and you
need just to define extra inputs for them. Then this will be exactly
what you have said, the Dirichlet will drive the solution. This what
ANSYS does and this is how I would solve this problem.

[tyler davis]

thank you for the suggestion

I want to update on my progress:

1. I have successfully set a non-zero initial state vector
(corresponding to 25°C)

2. I think there may be something wrong with my reduced state equation
matrices. Using the same MATLAB post-processing code, compare my
results:
http://picasaweb.google.com/lh/photo/Hy5cxCbxqPSjRn9quKBPIQ?feat=directlink

to those using sample files posted on your website (the Thruster
model): http://picasaweb.google.com/lh/photo/5ff5vS0GHEPnqwfTe44m5Q?feat=directlink

You can see that my model "blows up" with time, with all temperature
going towards infinity with increasing time. I noticed that the
thruster files take much longer to compute an equal time domain. Also
the A and E matrices for the thruster files are quite sparse, while
mine are completely filled. Can you think of a reason why my model
would go unstable?

3. Regarding the Dirichlet BCs - it appears necessary to split the K/A
matrix (conductivity) into Dirichlet nodes and non-Dirichlet nodes.
Does this need to happen prior to Model Reduction? I can't
intelligently modify the A matrix after it has been reduced since it
does not represent real physical nodes, only a generic state-space
variable

thank you,
Tyler

On Apr 2, 1:06 am, Evgenii Rudnyi <use...@rudnyi.ru> wrote:
> >> As for the not zero Dirichlet, this should be considered as an extra
> >> input. You first make Dirichlet zero and apply power dissipation. This
> >> will give you the first input. Then you make power dissipation zero and
> >> apply not zero Dirichlet equal to one. This gives you the second input
> >> that you can scale afterwards in the reduced model.
>
> >> It does not make sense to define Dirichlets as outputs, as they are known.
>
> > I agree that Dirichlet nodes are not outputs to the reduced model, but
> > rather should be inputs.  The temperature of the Dirichlet nodes
> > affects all the other temperatures in the reduced model, since they
> > are linked through conductors.  I don't understand your method, but I
> > don't think it will work for my model.  I want to use the reduced
> > model to evaluate different settings of the Dirichlet nodes (that
> > represent different environments).  I literally want to "plug in" the
> > defined temperatures of the boundary nodes, and have them propagate
> > through the other non-boundary nodes in the reduced model.
>
> > Perhaps it would help to show you a sample of the results from the
> > full-fidelity model (click on the link for image)

> >http://picasaweb.google.com/lh/photo/N4c8gDdQAj7NsisvbyjbxQ?feat=dire...

[Evgenii Rudnyi]

> 1. I have successfully set a non-zero initial state vector
> (corresponding to 25�C)
>
> 2. I think there may be something wrong with my reduced state equation
> matrices. Using the same MATLAB post-processing code, compare my
> results:
> http://picasaweb.google.com/lh/photo/Hy5cxCbxqPSjRn9quKBPIQ?feat=directlink
>
> to those using sample files posted on your website (the Thruster
> model): http://picasaweb.google.com/lh/photo/5ff5vS0GHEPnqwfTe44m5Q?feat=directlink
>
> You can see that my model "blows up" with time, with all temperature
> going towards infinity with increasing time. I noticed that the
> thruster files take much longer to compute an equal time domain. Also
> the A and E matrices for the thruster files are quite sparse, while
> mine are completely filled. Can you think of a reason why my model
> would go unstable?

The boundary conditions in the benchmark correspond to a very small film
coefficient. The stationary state does exist but it very high and not
physical. What you observe is presumably correct but here the sense
makes only simulation for about 1 s.

You will find the transient simulation from ANSYS to compare with at

http://evgenii.rudnyi.ru/doc/misc/pyros/thruster.tar.gz

see README there.

The matrices after the FEM discretization are sparse indeed. The
matrices in the reduced system however are dense.

> 3. Regarding the Dirichlet BCs - it appears necessary to split the K/A
> matrix (conductivity) into Dirichlet nodes and non-Dirichlet nodes.
> Does this need to happen prior to Model Reduction? I can't
> intelligently modify the A matrix after it has been reduced since it
> does not represent real physical nodes, only a generic state-space
> variable

Yes, this should be done before model reduction. MOR is applied to a
system of ordinary differential equations and the equations for
Dirichlets do not belong to them.

[tyler davis]

So do you think my inclusion of Dirichlet nodes (zero capacitance) in
my E matrix leads to the unstable behavior of the reduced model?

Are you recommending I remove the Dirichlet nodes from my A and E
system matrices, and apply their effects as heat loads (B matrix)?

Regarding scaling these loads after model reduction: the amount of
heat flowing between a Dirichlet node and a non-Dirichlet node is
dependent upon the current temperature of the non-Dirichlet node. So
the load magnitude cannot simply be scaled prior to a model run; it
load must somehow be updated during model simulation to account for
changing temperatures.

Do you have sample files that show thermal Dirichlet boundary
conditions being applied via the load vector (the way ANSYS handles
these loads)?

On Apr 3, 2:55 am, Evgenii Rudnyi <use...@rudnyi.ru> wrote:
> > 1. I have successfully set a non-zero initial state vector
> > (corresponding to 25 C)
>
> > 2. I think there may be something wrong with my reduced state equation
> > matrices.  Using the same MATLAB post-processing code, compare my
> > results:

> >http://picasaweb.google.com/lh/photo/Hy5cxCbxqPSjRn9quKBPIQ?feat=dire...

>
> > to those using sample files posted on your website (the Thruster

> > model):http://picasaweb.google.com/lh/photo/5ff5vS0GHEPnqwfTe44m5Q?feat=dire...

[Evgenii Rudnyi]

> So do you think my inclusion of Dirichlet nodes (zero capacitance) in
> my E matrix leads to the unstable behavior of the reduced model?

I do not know. I have just written that in the microtruster benchmark
the stationary state corresponds to very high temperature, say
temperature of the Sun.

> Are you recommending I remove the Dirichlet nodes from my A and E
> system matrices, and apply their effects as heat loads (B matrix)?

This would be more numerically robust. Well, I have to remember that
trick. I believe that long time ago I have tried MOR with it but I do
not remember now what the results were.

> Regarding scaling these loads after model reduction: the amount of
> heat flowing between a Dirichlet node and a non-Dirichlet node is
> dependent upon the current temperature of the non-Dirichlet node. So
> the load magnitude cannot simply be scaled prior to a model run; it
> load must somehow be updated during model simulation to account for
> changing temperatures.
>
> Do you have sample files that show thermal Dirichlet boundary
> conditions being applied via the load vector (the way ANSYS handles
> these loads)?

I guess I can find some. In what form do you want to obtain this? Is
ANSYS script okay?