Time Series Data Analysis: Logical segmentation of phases of data

I have high speed mixing data of various materials and during this mixing phase, the materials go through distinct phase changes. The super fancy mixer I have controls the rotor speed and trends the power required to achieve that speed.

During the last phase, you can see from the data attached, there appears to be 3 distinct phases in the required rotor power. The first phase, the material starts to clump, so the rotor requires a little more power, and as the material breaks up, we see a distinct decrease in power. As this is happening, the temperature is going up and at some point, the material densifies and granulates. When this happens, we see a sharp increase in the rotor power. This is very consistent but depending on what else I add to the mix, I can get different slopes, as some things delay granulation and other things speed them up.

So basically I have 100s of experiments that show these three distinct phases but all have very different time lengths. What I am trying to do is fit these three phases with very simple linear or polynomial fits so I can put more basic numbers into a design of experiment.

Since I haven't had a lot of time to look around I thought I would just throw this out there.

Attached is the data and some code that simply creates X linear segments and I am hoping some of you have some ideas on how I could create more logical code that would be able to find these three segments, for lack of a better term, that would then allow me to extract some more basic information.

Data and code attached. Sorry I haven't added any comments but you should be able to run this and see what I am trying to do.... Thanks

clc; clear all; 

load('EL1_Data.mat')
tic
% dSeg is how many brackets we are going to make
dSeg = 6;

% Create segments to linearly fit
brack = linspace(min(timeLast),max(timeLast), dSeg);
% Cell to store the different brackets
iX = cell(dSeg-1);
dataFits = zeros(dSeg-1,2);
%errors = struct();
figure(1)
plot(timeLast,rPowLast,'.')
hold on

for i = 1:dSeg-1
    iX{i} = find(timeLast >= brack(i) & timeLast <= brack(i+1));
    [p,s] = polyfit(timeLast(iX{i}),rPowLast(iX{i}),1);
    dataFits(i,:) = p;
    errors(i) = s;
    [yFit, yErr] = polyval(p,timeLast(iX{i}),s);
    plot(timeLast(iX{i}),yFit,'-')
    hold on

end
hold off
xlabel('Time (s)')
ylabel('Rotor Power (W)')

toc

[dpb]

Look at "piecewise linear regression". I demonstrated/provided code for the case of two segments not very long ago here (but, one of the problems Answers hasn't solved but was supposed to) I can't seem to find it easily again at the moment.

I've not, at least recently, moved beyond the two segments but it's just additional conditions on the model.