Support for partial derivatives of functions

[Adrian Pop]

It seems we have no support for partial derivatives of functions.
See the attached model and Modelica Spec 3.2 revision 2, Section 12.7.2.

The question is how we should handle this?
We save the information for the back-end?

We have a function definition such as:

function df_dqflow = der(flowCharacteristic, q_flow);

Use this script:

loadFile("tfm.mo"); getErrorString();
simulate(ThermoPower.Test.GasComponents.TestFanMech); getErrorString();

[Adrian Pop]

When we try to simulate this model we get:

ThermoPower/Gas.mo:3223:7-3226:21:writable] Error: Class df_dqflow not found in scope ThermoPower.Gas.FanMech (looking for a function or record).
Error: Error occurred while flattening model ThermoPower.Test.GasComponents.TestFanMech

[Adrian Pop]

Any comments/suggestions on how to handle partial derivatives of functions?
The specification seems to be rather vague about this.

[Francesco Casella]

My understanding is that Automatic Differentiation techniques should be used in order to process the algorithm of the function, resulting in the algorithm that computes the derivatives. Dymola has been doing that for 10 years now, we've used this feature in the past for terrain representation in vehicle simulation.

Please note I will remove these partial derivatives from ThermoPower right away (they are used for simplified models in homotopy, but I can figure out other simplified models that don't need this information). As long as you keep using Adrian's provided .mo file for testing, it's ok.

[Willi Braun]

Basically our differentiation model is also able to differentiate that, but we just haven't considered such cases so far.

Adrian, we need to call the differentiate module at a certain point in the frontend to generate a new differentiated function. So I would create an interface function therefore in the differentiation module.

[Adrian Pop]

Willi introduced a new interface Differentiate.differentiateFunctionPartial in r18683.
I'll use this in the front-end to generate a new function with the given name and
add it to the environment (a bit of a question on how here) and to the function tree.

[Martin Sjölund]

This ticket was not closed for 1.9.1, which has now been released. It was batch modified for milestone 1.9.2 (but maybe an empty milestone was more appropriate; feel free to change it).

[Martin Sjölund]

Milestone changed to 1.9.3 since 1.9.2 was released.

[Martin Sjölund]

Adrian: How did this work go?

[Adrian Pop]

No idea. Maybe Willi knows more.

[Adrian Pop]

Ah, we dö not yet handle it in the frontend but i guess should be easy to add.
I will have a look.

[Martin Sjölund]

Moved to new milestone 1.9.4

[Martin Sjölund]

Milestone pushed to 1.9.5

[Martin Sjölund]

Milestone renamed

[Martin Sjölund]

Ticket retargeted after milestone closed

[Martin Sjölund]

Milestone moved to 1.12.0 due to 1.11.0 already being released.

[Francesco Casella]

The milestone of this ticket has been reassigned to "Future".

If you think the issue is still valid and relevant for you, please select milestone 1.13.0 for back-end, code generation and run-time issues, or 2.0.0 for front-end issues.

If you are aware that the problem is no longer present, please select the milestone corresponding to the version of OMC you used to check that, and set the status to "worksforme".

In both cases, a short informative comment would be welcome.

[Francesco Casella]

Simpler test model:

model TestDerivative
  function f
    input Real x;
    input Real y;
    output Real z;
  algorithm
    z := 1;
    z := z + x*y;
  end f;
  
  function dfdx = der(f,x);
  
  Real x = time;
  Real df;
equation
  df = dfdx(x, x^2);
end TestDerivative;

Gives the following error with both OF and NF:

[1] 01:14:09 Translation Error
[TestDerivative: 16:3-16:20]: Class dfdx not found in scope TestDerivative
(looking for a function or record).

[Francesco Casella]

I think this should be resolved in 2.0.0. We do have all the needed primitives in the backend, as I understand it's just a matter of connecting the dots from the NF.

[Francesco Casella]

@perost, can you coordinate with @Karim and @AnHeuermann to get this done?

It should be rather straightforward.

[Per Östlund]

Replying to casella:

It should be rather straightforward.

That depends on what kind of support we need for partial derivatives of functions in the NF. The Differentiate module is made for the old DAE structures, so using it in the NF is not really possible (the NF can translate to OF format, but too much is lost in that translation which makes translating back to NF format very nontrivial).

This means that if we want the NF to be able to generate the body of a partial derivative we'd have to rewrite the Differentiate module using the NF structures. This is needed if we want the NF to be able to evaluate partial derivatives of functions.

If we don't care about being able to evaluate such functions I guess we can just pass them through the frontend, in which case we only need to know their type (which is the same as the function being derived?).

Then we can generate the partial derivatives either during the conversion to DAE or somewhere in the backend. Differentiate.differentiatePartialFunction takes some backend structures as argument that the NF can't generate, so I guess where to do it depends on how important those arguments are.

[Francesco Casella]

Replying to perost:

If we don't care about being able to evaluate such functions

I think this is definitely the case. Of course it may be good for efficiency reason to be able do so in some cases, but if we can't, it's no big deal.

I guess we can just pass them through the frontend,

That was my idea. The machinery to compute total derivatives of functions is only there, computing partial derivatives is actually a special case of that, so it should be easy to implement.

in which case we only need to know their type (which is the same as the function being derived?).

Not exactly, the unit of the differentiated funtion output is of course the unit of the original function output divided by the unit of the variable w.r.t. which the function is differentiated. As far as I understand, everything else should be the same.

Then we can generate the partial derivatives either during the conversion to DAE or somewhere in the backend. Differentiate.differentiatePartialFunction takes some backend structures as argument that the NF can't generate, so I guess where to do it depends on how important those arguments are.

Never mind that, let's leave the job to the backend, you only need to prepare the patient on the operating table.

[Karim Adbdelhak]

Differentiating a function w.r.t. a cref should not be a problem, as long as we know all the dependencies. But i don't really understand what is supposed to happen in your example.

If we differentiate the function

  function f
    input Real x;
    input Real y;
    output Real z;
  algorithm
    z := 1;
    z := z + x*y;
  end f;

like you imply with this call w.r.t x

  function dfdx = der(f,x);

It would look somewhat like this:

  function dfdx
    input Real x;
    input Real y;
    output Real z;
  algorithm
    z := 0;
    z := z + y;
  end dfdx;

so the equation

  df = dfdx(x, x^2);

could be simplified to be

  df = x^2;

Is that what you would expect?

EDIT: Scrap that! I forgot the actual partial derivative of the input:

  function dfdx
    input Real x;
    input Real y;
    output Real z;
  algorithm
    z := 0;
    z := z + y + x * der(y,x);
  end dfdx;

I guess i have to put some thought into this... what do you expect happening here? :D

[Francesco Casella]

I'm not an expert in AD, but what I understand is that you keep the original algorithm and add the instructions to compute the derivatives alongside it. Hence, if you take

  function f
    input Real x;
    input Real y;
    output Real z;
  algorithm
    z := 1;
    z := z*x + z*y + x*y;
    z := z*z;
  end f;

then der(f,x) can be computed as

  function df
    input Real x;
    input Real y;
    output Real dz_x;
  algorithm
    dz_x := 0;
    z    := 1;

    dz_x := z   +   0 + y
    z    := z*x + z*y + x*y;

    dz_x := 2*z*dz_dx;
    z    := z*z;
  end f;

Indeed, if you simbolically substitute all the steps, the first function returns (x + y + x*y)^2, whose derivative w.r.t. x is 2*(x + y + x*y)*(1+y).

Then, for example, df_dx(x = 2, y = 3) = 2*(2 + 3 + 6)*(1 + 3) = 88, which is the same result you get when running the differentiated function:

dz_x := 0;
z := 1;
dz_x := 1 + 3 := 4;
z := 1*2 + 1*3 + 2*3 = 11;
dz_dx := 2*11*4 := 88
z := 2*11 = 22; // this is actually useless

Mind you, there is no der(x,y) term, those are two independent variables of the function.

I'm also not expecting that the partial derivative function is further simplified or inlined or symbolically processed in any way. It should simply be compiled into executable code that returns the appropriate values.

[Francesco Casella]

​https://en.wikipedia.org/wiki/Automatic_differentiation

[Karim Adbdelhak]

Yea you are totally right!

[Francesco Casella]

1.16.0 feature freeze is imminent, new developments postponed to 1.17.0

[Francesco Casella]

Replying to perost:

If we don't care about being able to evaluate such functions I guess we can just pass them through the frontend, in which case we only need to know their type

Yes! There will never be the need of evaluating such functions in the front end in practice.

(which is the same as the function being derived?).

Yes, except that the unit of the output is the unit of the original output divided by the unit of the variable you are differentiating with respect to

[Francesco Casella]

@Karim, I will need this feature for the project dealing with the CO2 mixture and Peng-Robinson equation, do you think you can do that?

[Francesco Casella]

Rescheduled to 1.18.0

[Francesco Casella]

Ticket retargeted after milestone closed

[Francesco Casella]

1.18.0 blocker tickets moved to 1.19.0