Embedded Matlab in Simulink to find a zero by bisection method
Red Star
in a larger Simulink program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function b=find_b(Gamma_0,k,p_inf)
% This block supports the Embedded MATLAB subset.
% See the help menu for details.
b1=0.01;
b2=10;
while abs(b1-b2)>0.01
f1=-1/4*b1+k*b1^2/4*p_inf-k*b1^2/4+k*b1/2+b1/2-(Gamma_0^2./
(8*pi^2.*b1))/k;
f2=-1/4*b2+k*b2^2/4*p_inf-k*b2^2/4+k*b2/2+b2/2-(Gamma_0^2./
(8*pi^2.*b2))/k;
bc=(b1+b2)/2;
f=-1/4*bc+k*bc^2/4*p_inf-k*bc^2/4+k*bc/2+bc/2-(Gamma_0^2./
(8*pi^2.*bc))/k;
if f>0 && f1>0
b1=bc;
elseif f>0 && f2>0
b2=bc;
end
if f<0 && f1<0
b1=bc;
elseif f<0 && f2<0
b2=bc;
end
end
b=bc;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ok...maybe my script doesn't do what I would, anyway...why Matlab
isn't able to assign a value to b and stop my simulation? bc is
defined in 'while' cycle!
Maybe 'while' cycle doesn't work in Embedded enviroment.
If I write a simple M-file, out of Simulink, I get an output (not the
correct result because my script is wrong, but it runs! )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Anyone know how to write a Embedded Matlab subset to find a zero by
bisection method?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
someone
% Not exactly sure, but shouldn't b1 & b2 be swapped in second if-else-end condition?
if f>0 && f1>0
b1=bc;
elseif f>0 && f2>0
b2=bc;
end
if f<0 && f1<0
b1=bc;
elseif f<0 && f2<0
b2=bc;
end
end
% should be:
if f>0 && f1>0
b1=bc;
elseif f>0 && f2>0
b2=bc;
end
if f<0 && f1<0
b2=bc;
elseif f<0 && f2<0
b1=bc;
end
end
% or, perhaps even better:
if f>0
if f1>0
b1=bc;
elseif f2>0
b2=bc;
end
elseif f<0
if f1<0
b2=bc;
elseif f2<0
b1=bc;
end
end
% What should happen if f or f1 == 0?
% Why not use MATLAB fzero function?
Michael Hosea
pre-define b before the loop because the compiler doesn't notice that the
loop will have to execute at least once, but once you get the algorithm
working correctly in MATLAB *and* get a few little things squared away like
defining b in front of the loop, it will work.
--
Mike
"Red Star" <matteo.d...@gmail.com> wrote in message
news:54a60594-f284-4244...@v15g2000hsa.googlegroups.com...
Red Star
% What should happen if f or f1 == 0?
% Why not use MATLAB fzero function?
I thought it's impossible f or f1 == 0 exactly.
I cannot use 'fzero' because Embedded Matlab subset doesn't accept
anonymous function (and fzero needs that)
Red Star
Mike, I didn't understand where is my error.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function b = find_b(k,p_inf,Gamma_0,b1,b2)
% This block supports the Embedded MATLAB subset.
% See the help menu for details.
bc=0;
while abs(b1-b2)>0.1
f1=-1/4*b1+k*b1^2/4*p_inf-k*b1^2/4+k*b1/2+b1/2-(Gamma_0^2./
(8*pi^2.*b1))/k;
f2=-1/4*b2+k*b2^2/4*p_inf-k*b2^2/4+k*b2/2+b2/2-(Gamma_0^2./
(8*pi^2.*b2))/k;
bc=(b1+b2)/2;
f=-1/4*bc+k*bc^2/4*p_inf-k*bc^2/4+k*bc/2+bc/2-(Gamma_0^2./
(8*pi^2.*bc))/k;
if f>0 && f1>0
b1=bc;
elseif f>0 && f2>0
b2=bc;
end
if f<0 && f1<0
b1=bc;
elseif f<0 && f2<0
b2=bc;
end
end
b=bc;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Now a 'b' value is given but it's simply the difference between b1 and
b2
b1=50
b2=0.1
The while condition, 'abs(b1-b2)<0.1' is satisfied, but the script
doesn't seem go into the cycle.
I'd like a infinite loop than this dummy result !
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% How would you correct my script to find the zero of the 'function
f' in [b1 b2] interval by a Embedded Matlab subset? %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Michael Hosea
algorithm. You never said what values of k, p_inf, and Gamma_0 you were
using.
--
Mike
"Red Star" <matteo.d...@gmail.com> wrote in message
news:193fc800-4a1f-4ebc...@b1g2000hsg.googlegroups.com...
Red Star
> I would write it quite differently, but I don't think your problem is the
> algorithm. You never said what values of k, p_inf, and Gamma_0 you were
> using.
These values come from the simulation (they are inputs of the block).
I'm not able to foresee them, but I'm quite sure they are in a range
that don't give me problems.
The only trouble is that these parameters are variable in the time,
but I would that b is freezing at, let's say, 't=t*-eps' when t>t*.
That is in the initial part of my simulation the script must find 'b'
and it has to hold that until the simulation's end (I wrote a post
regarding).
This aspect is more important: I could cheat and simplfy 'f', so I'd
determine 'b' in a close form. Anyaway finding the correct 'b' is to
prefer.
Do you suggest a solider algorithm? Maybe a Newton's method is better?
Or else there is a pre-built algorithm in Matlab? (fzero seems to not
work in Embedded version)
Michael Hosea
worked. I can't tell whether there is any issue with Embedded MATLAB (e.g.
some bug in an earlier version) if I don't have all the parameters for at
least *one* test case where you think it is giving you the wrong answer.
Bisection is robust but not the fastest algorithm. Given the low accuracy
you are requiring here, I think it is fine, provided you don't have any
difficulty providing the initial interval bracketing a sign change in f(x).
That might be a good thing to check if you have any doubts at all that f(b1)
and f(b2) have opposite signs.
It is possible to call fzero, but you have to put the call site and the
function definition in an external M file that is called from the Embedded
MATLAB block and declared extrinsic. The drawback, aside from needing to
understand what I just wrote, is that you couldn't generate code from that
using RealTime Workshop.
--
Mike
"Red Star" <matteo.d...@gmail.com> wrote in message
news:e2598bf7-2555-436f...@25g2000hsx.googlegroups.com...
Red Star
I wrote a simple script to test 'f'
function [bc,f]=testb(k,p_inf,Gamma_0,b1,b2)
bc=b1:0.1:b2;
f=-1/4*bc+k.*bc.^2/4*p_inf-k.*bc.^2/4+k*bc/2+bc/2-(Gamma_0^2./
(8*pi^2.*bc))/k;
plot(bc,f);
grid on
end
with:
testb(4,1,0.825,-200,1000)
Zero obviously is a singularity point.
for bc<0 we have f<0 and for bc>0 we have f>0.
This behavior is more evident using: [x y]=testb(4,1400,1000,-10,10);
So the problem is that there isn't any root of 'f'. Maybe there is a
combination fo parameters that give a root a 'f', I don't know...
anyway I suppose I have to find another method to estimate 'b'.
(probably I should do other assumptions in my model)
The only trouble persisting is about 'freezing' my result.
My block has k, p_inf, Gamma_o like inputs.
'b' is the only output.
k, p_inf are constants.
Gamma_0 is a time function.
Hence 'b' is a time function. The block calculates a different b at
each t, of course.
At t=t* I would 'b' doesn't change anymore, but it stays constant with
the last calculated value until the end of my simulation.
I'd have this behavior:
let's suppose t*=10 and simulation lasts 18
t Gamma_0 b
----------------------------------
3 0.8 1
6 0.9 0.3
9 0.7 1.5
12 1.2 1.5
15 0.5 1.5
18 0.9 1.5
(The relation between 'Gamma_0' and 'b' is not importnat in this
example)
I tried to use 'switch' and 'workspace', but values stored in
workspace are available only at end of the simulation.
There is a pre-built block solving my problem, or how can I do ?
Thank you Mike
Matteo
someone
All I was suggesting is that one of the if-else-end conditions should test for f >=0 (or f <= 0) to make your code more robust. (Or maybe you should terminate the loop if f == 0.) How do you know f == 0 is impossible? f==0 may be an event that occurs with probability zero, but its not necessarily impossible.
The function input to fzero can be an m-file fuction.
Red Star
> Red Star <matteo.diplom...@gmail.com> wrote in message <fae8e57c-70e6-45f6-8335-f3515cd0f...@k7g2000hsd.googlegroups.com>...
Yes, you're right.
But what can you tell me about my second trouble? (see the example in
my last thread)
someone
> 6-8335-f3515cd0f...@k7g2000hsd.googlegroups.com>...
> >
> > > % What should happen if f or f1 =3D=3D 0?
> > > % Why not use MATLAB fzero function?
> >
> > > I cannot use 'fzero' because Embedded Matlab subset doesn't accept
> > > anonymous function (and fzero needs that)
> >
> t for f >=3D0 (or f <=3D 0) to make your code more robust. =A0(Or maybe you=
> should terminate the loop if f =3D=3D 0.) =A0How do you know f =3D=3D 0 is=
> impossible? =A0f=3D=3D0 may be an event that occurs with probability zero,=
> but its not necessarily impossible.
> >
> > The function input to fzero can be an m-file fuction.
>
> Yes, you're right.
> But what can you tell me about my second trouble? (see the example in
> my last thread)
If I understand it, your second trouble is that b (bc inside the while loop) is not changing after some number of iterations. One way this can happen is if the while loop intervals (f or f1 delta values) are <= eps. (Due to machine round-off problems. Then the bisection method simply oscillates between +/- eps.) In this case, terminate the loop when bc does not change from its previous value. (This is as close to zero as you can get.)
Or did I misinterpret your second trouble?
The MATLAB fzero function is simply an m-file. Call it up in the editor and see how MathWorks handles this condition. Last time I looked, its implementation of the bisection method was pretty similar to yours.
Red Star
> Red Star <matteo.diplom...@gmail.com> wrote in message <d1faa1c5-8670-4be6-8c1b-898894599...@t41g2000hsc.googlegroups.com>...
> > On 3 Ott, 17:17, "someone " <some...@somewhere.net> wrote:
> > > Red Star <matteo.diplom...@gmail.com> wrote in message <fae8e57c-70e6-45f=
> > 6-8335-f3515cd0f...@k7g2000hsd.googlegroups.com>...
>
> > > > % What should happen if f or f1 =3D=3D 0?
> > > > % Why not use MATLAB fzero function?
>
> > > > I thought it's impossible f or f1 =3D=3D 0 exactly.
> > > > I cannot use 'fzero' because Embedded Matlab subset doesn't accept
> > > > anonymous function (and fzero needs that)
>
> > > All I was suggesting is that one of the if-else-end conditions should tes=
> > t for f >=3D0 (or f <=3D 0) to make your code more robust. =A0(Or maybe you=
> > should terminate the loop if f =3D=3D 0.) =A0How do you know f =3D=3D 0 is=
> > impossible? =A0f=3D=3D0 may be an event that occurs with probability zero,=
> > but its not necessarily impossible.
>
explained my problem too approximatly. However I solved that.
Anyway...you can post a *working* Embedded Matlab subset that find a
zero?
Many people suggest me to call an external M-file with fzero routine,
but I am not able to do.
My problem is: finding the root of a 'f' function, with some
parameters that can change whenever and they must be inputs of the
routine.
What do you suggest I could do,please?
someone
> 6-8c1b-898894599...@t41g2000hsc.googlegroups.com>...
> > > On 3 Ott, 17:17, "someone " <some...@somewhere.net> wrote:
> -45f=3D
> > > 6-8335-f3515cd0f...@k7g2000hsd.googlegroups.com>...
> >
> > > > > % What should happen if f or f1 =3D3D=3D3D 0?
> > > > > % Why not use MATLAB fzero function?
> >
> > > > > I cannot use 'fzero' because Embedded Matlab subset doesn't accept
> > > > > anonymous function (and fzero needs that)
> >
> tes=3D
> > > t for f >=3D3D0 (or f <=3D3D 0) to make your code more robust. =3DA0(Or=
> maybe you=3D
> > > =A0should terminate the loop if f =3D3D=3D3D 0.) =3DA0How do you know f=
> =3D3D=3D3D 0 is=3D
> > > =A0impossible? =3DA0f=3D3D=3D3D0 may be an event that occurs with proba=
> bility zero,=3D
> > > =A0but its not necessarily impossible.
> >
> The second trouble is not that. You didn't understarnd because I
> explained my problem too approximatly. However I solved that.
>
> Anyway...you can post a *working* Embedded Matlab subset that find a
> zero?
> Many people suggest me to call an external M-file with fzero routine,
> but I am not able to do.
>
> My problem is: finding the root of a 'f' function, with some
> parameters that can change whenever and they must be inputs of the
> routine.
> What do you suggest I could do,please?
% write a function m-file and save it as 'myfun.m'
function y = myfun(x)
y = x.^3-2*x-5;
% then use fzero to find the zero near 2
z = fzero('myfun',2)
z =
2.0946
% Now, if you use more than one input to the function:
function y = myfun(x,a)
y = a + x.^3-2*x-5;
% You will need to use feval
Red Star
> Red Star <matteo.diplom...@gmail.com> wrote in message <f0a788b1-cb79-4b27-93fd-f66fb89f3...@w7g2000hsa.googlegroups.com>...
Ok, this is the solution purposed by the manual.
But it works in Embedded mode?Did you try?And what means 'I will need
to use feval'?
I need a example to understand.
Let's suppose my function is 'f'
I should writhe the f.m :
function y=f(x,A,B,)
y=x.^2+A*x+B;
end
so in simulink I need an Embedded Matlab subset calling this f.m. How
must I do?
I'd create a block (the Embedded Matlab subet) with inputs A B and y0
(root's 'f') as output.
But what must I write in that script, inside the block?
Let's remember that is an Embedded Matlab subset and not a M-file:
some differences could be.
Michael Hosea
algorithm within an Embedded MATLAB block. However, if you don't need to
generate code with RealTime Workshop, here is how to call a "FunFun" like
FZERO from an Embedded MATLAB block:
%%% Begin Embedded MATLAB Block.
function y = fcn(a,b) %#eml
% This function will work in Simulink, but it will not work with RealTime
% Workshop code generation.
eml.extrinsic('foo');
y = 0; %#ok<NASGU>
y = foo(1,a,b); % Call MATLAB
%%% End Embedded MATLAB Block.
%%% Begin foo.m, an ordinary MATLAB M file to reside
%%% on the MATLAB path outside of model.
function r = foo(guess,a,b)
r = fzero(@(x)f(x,a,b),guess);
function y = f(x,a,b)
y = sin(a*x + b); % or whatever function you want.
%%% End extrinsic M file.
Note that Embedded MATLAB supports function handles, but not *anonymous*
function handles. You can pass function handles around within an Embedded
MATLAB block, but you cannot pass them or otherwise refer to functions
defined in Embedded MATLAB to any extrinsic function. This is why I wrote
the "wrapper function" foo.m, which just accepts the *parameters* required
to define the function. An extrinsic function, by the way, is one that
that runs in the MATLAB environment just as it would when given the same
inputs from the MATLAB command prompt. Consequently, any solution requiring
an extrinsic function cannot be used for Realtime Workshop code generation.
--
Mike
"Red Star" <matteo.d...@gmail.com> wrote in message
news:c972032b-bb02-487c...@z66g2000hsc.googlegroups.com...