Is there a way to improve the code speed?

Hi,

Here below is the code I made for plotting the output frequency over time of a classic Phase Lock Loop module.

The code is horrible in speed, much slower than the simulink model which I also attached here.

My team leader coded similar thing in Julia but in a quite different way, while the speed is ultra fast. Within less than 3 seconds can he get the plot. I have realized there is something like VectorContinuousCallback in the Julia code for solving ODEs and many other special things.

Can we improve the code speed in matlab although without such functions?

clear;close all;clc;
tr(1)=5e-6;tref=25e-9;tdd=250e-12;ip=60e-6;in=60e-6;c0=40e-12;c1=360e-12;c2=20e-12;r1=14000;r2=1000;f0=14.16e9;kvco=300e6;n=360;
m=2500;ic=zeros(1,m);tc=sym('tc',[1 m]);p=sym('p',[1 m]);q=sym('q',[1 m]);s=sym('s',[1 m]);phy=sym('phy',[1 m]);syms f(t) g(t) h(t) vc0(t) vc1(t) vc2(t) t x a b;
ode=[r1*c1*diff(vc1,t)+vc1==vc0, r2*c2*diff(vc2,t)+vc2==vc0, c0*diff(vc0,t)+c1*diff(vc1,t)+c2*diff(vc2,t)==0];conds=[vc0(0)==0,vc1(0)==0,vc2(0)==0];f(t)=dsolve(ode,conds).vc0;g(t)=dsolve(ode,conds).vc1;h(t)=dsolve(ode,conds).vc2;
t0=vpasolve(mod(0,1)+int((f0+kvco*h)/n,0,x)==1,x,[0,inf]);
if t0==tr(1)
    ic(1)=ip-in;tc(1)=tr(1);p(1)=f(tc(1));q(1)=g(tc(1));s(1)=h(tc(1));phy(1)=0;
elseif t0>tr(1)
    ic(1)=ip;tc(1)=tr(1);p(1)=f(tc(1));q(1)=g(tc(1));s(1)=h(tc(1));phy(1)=int((f0+kvco*h)/n,0,tc(1));
elseif t0<tr(1)
    ic(1)=-in;tc(1)=t0;p(1)=f(tc(1));q(1)=g(tc(1));s(1)=h(tc(1));phy(1)=0;
end
fplot(f0+kvco*h,[0,double(tc(1))]);hold on
for i=2:m
    if ic(i-1)==-in
        ic(i)=ip-in;
        if tc(i-1)-tr(1)>=0
            tc(i)=tr(1)+tref*(floor((tc(i-1)-tr(1))/tref)+1);
        else
            tc(i)=tr(1);
        end
        ode=[r1*c1*diff(vc1,t)+vc1==vc0, r2*c2*diff(vc2,t)+vc2==vc0, c0*diff(vc0,t)+c1*diff(vc1,t)+c2*diff(vc2,t)==ic(i-1)];conds=[vc0(0)==p(i-1),vc1(0)==q(i-1),vc2(0)==s(i-1)];f(t)=dsolve(ode,conds).vc0;g(t)=dsolve(ode,conds).vc1;h(t)=dsolve(ode,conds).vc2;
        p(i)=f(tc(i)-tc(i-1));q(i)=g(tc(i)-tc(i-1));s(i)=h(tc(i)-tc(i-1));phy(i)=phy(i-1)+vpaintegral((f0+kvco*h)/n,0,tc(i)-tc(i-1));
    elseif ic(i-1)==ip
        ic(i)=ip-in;
        ode=[r1*c1*diff(vc1,t)+vc1==vc0, r2*c2*diff(vc2,t)+vc2==vc0, c0*diff(vc0,t)+c1*diff(vc1,t)+c2*diff(vc2,t)==ic(i-1)];conds=[vc0(0)==p(i-1),vc1(0)==q(i-1),vc2(0)==s(i-1)];f(t)=dsolve(ode,conds).vc0;g(t)=dsolve(ode,conds).vc1;h(t)=dsolve(ode,conds).vc2;
        tc(i)=tc(i-1)+vpasolve(mod(phy(i-1),1)+int((f0+kvco*h)/n,0,x)==1,x,[0,inf]);
        p(i)=f(tc(i)-tc(i-1));q(i)=g(tc(i)-tc(i-1));s(i)=h(tc(i)-tc(i-1));phy(i)=0;
    elseif ic(i-1)==ip-in
        if tc(i-1)-tr(1)>=0
            a=tr(1)+tref*(floor((tc(i-1)-tr(1))/tref)+1);
        else
            a=tr(1);
        end
        ode=[r1*c1*diff(vc1,t)+vc1==vc0, r2*c2*diff(vc2,t)+vc2==vc0, c0*diff(vc0,t)+c1*diff(vc1,t)+c2*diff(vc2,t)==ic(i-1)];conds=[vc0(0)==p(i-1),vc1(0)==q(i-1),vc2(0)==s(i-1)];f(t)=dsolve(ode,conds).vc0;g(t)=dsolve(ode,conds).vc1;h(t)=dsolve(ode,conds).vc2;
        b=tc(i-1)+vpasolve(mod(phy(i-1),1)+int((f0+kvco*h)/n,0,x)==1,x,[0,inf]);
        if b==a
            ic(i)=ip-in;tc(i)=a;p(i)=f(tc(i)-tc(i-1));q(i)=g(tc(i)-tc(i-1));s(i)=h(tc(i)-tc(i-1));phy(i)=0;
        elseif b>a
            ic(i)=ip;tc(i)=a;p(i)=f(tc(i)-tc(i-1));q(i)=g(tc(i)-tc(i-1));s(i)=h(tc(i)-tc(i-1));phy(i)=phy(i-1)+vpaintegral((f0+kvco*h)/n,0,tc(i)-tc(i-1));
        else
            ic(i)=-in;tc(i)=b;p(i)=f(tc(i)-tc(i-1));q(i)=g(tc(i)-tc(i-1));s(i)=h(tc(i)-tc(i-1));phy(i)=0;
        end
    end
    fplot(f0+kvco*h(t-tc(i-1)),[double(tc(i-1)),double(tc(i))])
end

[Kai Wang]

Matlab has no such function as WhenEvent is the ode solver in Mathematica or VectorContinousCallback in Julia.

We have to define the events outside the solver and run cycle to cycle.

With WhenEvent, the speed can be reduced to 3 minutes in Mathematica, very huge improvement.

[Dyuman Joshi]

It's not clear to me what you are trying to do and it's hard to follow the code. So I will ask some questions -

Any particular reason why you have defined tc, phy, p, q, and s as a symbolic variable and not numerical like ic?

The plots are added to the same figure, so they don't seem to be the bottle neck.

Have you tried using profiler on your code?

[Kai Wang]

Thank you for your questions.

I defined those parameters as symbolic variables because I assumed they would provide better accuracy in vpasolve. Maybe it is unnecessary.

I did not know such thing as profiler. Thank you for reminding, let me check with that.

Answer 1

[Sulaymon Eshkabilov]

Without digging much into your code, I'd suggest to put:

fplot(f0+kvco*h(t-tc(i-1)),[double(tc(i-1)),double(tc(i))]) outside of the loop that speeds up the simulation of your code signifcantly. Because you are plotting calculated data for over 2500 times <-- m=2500 which requires lots of time.

All the best.

[Kai Wang]

Thank you, that is a good point. I tried to comment the fplot, the speed merely changed.

The time consuming part is the ode solving.