How to compute numerical gradient of an unknown function in matlab Simulink?

Hi! I have a function that is unknown and that is being generated at every time instant t in my simulink model (Imagine it like a black box, wherein, I have data for at every time instant but not the mapping itself, the data looks something like this: ). In order to implement my algorithm, I need to find the direction that corresponds to greatest decrease (or increase) of this scalar mapping at every time instant. So, how can I find the numerical gradient descent of this unknown function at every time step, i.e at any ? The only thing known is that is lipshitz continuous. Kindly suggest me ways in which it can be implemented in simulink. Thanks for your time and consideration.

Answer 1

[Torsten]

The i'th component gi of the gradient vector g is approximately given by

gi = (f(x1,...,xi+h,...,x10)-f(x1,...,xi,...,x10))/h

Thus for each time t, given (x1,...,x10), you will have to evaluate f at (x1+h,x2,...,x10),(x1,x2+h,x3,...,x10),...,(x1,x2,...,x10+h) and calculate the gi.

[Amardeep Mishra]

The thing is, h will depend on the sampling time of the simulation. So if lets say I choose very very small sampling time, maybe something like , then in that case, I might try:

Will this term always exist?

I hope this is what you're suggesting me to do? Also as a sidenote, how can I use gradient in this case, because apparently someone suggested gradient for a similar problem here:

https://stackoverflow.com/questions/29478838/whats-the-best-way-to-calculate-a-numerical-derivative-in-matlab

[Torsten]

The sampling time in t is irrelevant for the gradient computation.

For each t, you have a vector

x(t)=(x1(t),...,x10(t)).

To compute the gradient of f in x(t), you can choose a value h>0 and evaluate

fh1 = f(x1(t)+h,x2(t),...,x10(t)), 
fh2 = f(x1(t),x2(t)+h,...,x10(t)),
...
fh10 = f(x1(t),x2(t),...,x10(t)+h).

With

f0 = f(x1(t),x2(t),...,x10(t)), 

the gradient gradf of f at x(t) is approximately

gradf = ((fh1-f0)/h,(fh2-f0)/h,...,(fh10-f0)/h)

[Amardeep Mishra]

I wish to point out that it is not possible to do the following because is unknown and is only generated at each time step.

fh1 = f(x1(t)+h,x2(t),...,x10(t)), 
fh2 = f(x1(t),x2(t)+h,...,x10(t)),
...
fh10 = f(x1(t),x2(t),...,x10(t)+h).

[Torsten]

But you get values for f(x1(ti),x2(ti),...,x10(ti)). How is f given ? By a differential equation ?

[Amardeep Mishra]

yes! is nothing but , where and ,

where :

(a) is the state of the system.

(b) is known as drift dynamics and

(c) is what we call as input coupling dynamics.

It is the control input that actually depends on w.

So, in short, yes the quantity Σ depends on both x and . In order to obtain x, one must integrate the differential equation: .

In other words, I need to find the direction that corresponds to maximum increase or decrease of the time series with respect to the parameter space

[Torsten]

And how do you get w and u(w) ? Are these explicit functions of t ?

[Amardeep Mishra]

w is computed based on a differential equation, that is there is an expression for wdot, upon whose integration one obtains, w and based on that u(x,w) is generated

[Torsten]

Then, since

f(w1,w2,...,w10) = x'*(p(x)+q(x)*u(x,w1,w2,...,w10)

I don't understand why you write that

fh1 = f(w1+h,w2,...,w10),

fh2 = f(w1,w2+h,...,w10),

...

fh10 = f(w1,w2,...,w10+h)

cannot be evaluated to compute the gradient of f with respect to w.

[Amardeep Mishra]

Because that would mean that, we'd have to perturb the parameters, and then perturb u and then integrate xdot to find the next state. Now at any given instant of time you can either give u(w). How can you give both u(w) and u(w+h) to the system?

[Torsten]

How can you give both u(w) and u(w+h) to the system?

If this is a technical problem that concerns SIMULINK, I can't answer it.

As you wrote, u is the control input specified by you. So in general, it should be possible to evaluate it at (x(t),w(t)) and (x(t),perturbed w(t)).