@vv Thanks for explanation

Hello everyone again,

I repeated the identification of the DC motor. In this case I used different DC servo motor and good measurements device.  Bas on collected data I obtained DC motor parameters in matlab. Below the results

First chart is the position of rotor, second speed and last is current of rotor . I modified model in simulink. I set the minimum and maxiumum value of integrator block in simulink becasue for exampel in secend chart the sposition is from 10 to -10V. All unitis in chart are volts. Leiter I will chcange this usnits to rad, rad/s and A.

Below csv file with data:

Column 1) - time

Column 2) - Input signal

Column 3) Position

Column 4) Speed

Column 5) current

DC_DATA.zip

@Preben Alsholm  Thank you. As I wrote above I will try to measure the i(t) in next week. 

@acer thanks. I will try to measure the i(t) in next week. This will help to find all parameters of dc motor.

@Preben Alsholm Thank you. We can not eliminate i(t). We need to estimate all parameters because at the end we change the representation of system to state space form 

A = [0      1       0
     0      -b/J    K/J
     0      -K/L    -R/L];

B = [0;...
     0;
     1/L];
C = [1   0  0;
     0   1  0];

And the typical value of above parameters are as I was show in my previous post:

J = 6.9347e-005;
b = 2.2480e-004;
K = 0.0094;
R = 0.2124;
L = 0.0015;

It means that: J, b should be small, L and K also and R can be bigger for example 1,2,4. We need to look on physic meaning of this parameters.

@acer  Thank you. Now I need a reasonable fit and I need a code/solution of my problem. As You see, results from matlab are not perfect. I need the results similar for results from matlab.

The position of dc motor was measured by incremental encoder so the chart is very smooth. The seed of dc motor was measured by tachogenerator with poor quality of sensor. This is the reason why we see a huge measurements noise.

Because at first I did not explain my problem exactly I decided to provide a detailed explanation below to clear the discussions.

1) I have dc motor model (differential equation) from here: http://ctms.engin.umich.edu/CTMS/index.php?example=MotorPosition&section=SimulinkModeling

2) I used this ODE eqution and model in simulink to find the prarmteres:

• J   - moment of inertia of the rotor 
• b  -  motor viscous friction constant
• K -  the motor torque
• R - electric resistance 
• L - electric inductance

To do this I measured the real speed of dc motor and position of rotor using PC computer and measurement card. How I messured this varibalies ?.I change the voltage applied to dc motor using measurement card of course. My impu signal V was a sinusoidal signal of parametres:

• Amplitude 0.5 volts
• Bias 0.5 volts

Samplig Time 0.01[s]. The input signal is as follows 

Ok, the respons of systems 

Ok, next I used matlab and simulink to find the parameters J, b, K, R, L of my model using data from csv file. This data was collected using measurement card. In this file I have data from measurement card

The columns in this file are:

• Column 1 - Is time vector (a read data with sampling time 0.01[s] for 0[s] to 50[s]
• Column 2 - Is input signal for dc motor V. THIS SIGNAL I GENERATE USING BLOK IN SIMULINK.
• Column 3 -Is position of rotor theta
• Column 4 - Is speed of rotor.

I did not measure the current i(t) because I can't measure this using measurement card. I need a different device. So now I only use wath I have to find parmatres of system. 

What next? Next I change my model in simulink as follows:

Finally I wrote the code to find the parameters of my model:

clc,close all,clear
%% Settings
PlotCSVData = 0;
UsePreviousParameters = 1;
%% Load data from csv file
LoadCSVData('servo_open_ident_1_A.csv')

%% Time
t = servo_open_ident_1_A(:,1);

%% Input
In(1).values = servo_open_ident_1_A(:,2);

%% Output
Out(1).values = servo_open_ident_1_A(:,3);
Out(2).values = servo_open_ident_1_A(:,4);

%% Plot open data
if PlotCSVData == 1
    figure;
    set(gcf,'Position',[291 277 1325 597])
    for i=1:length(In)
        subplot(length(In)+ length(Out),1,i)
            plot(t,In(i).values);
            xlabel('Time[s]');
            ylabel(['In(',num2str(i),')'])
            grid on
    end
    for i=length(In)+1:(length(In)+length(Out))
        subplot(length(In)+ length(Out),1,i)
            plot(t,Out(i-1).values);
            xlabel('Time[s]');
            ylabel(['Out(',num2str(i-1),')'])
            grid on
    end
end
%% Ident - Set defult value of paramteres and I/O data
if UsePreviousParameters == 0
    J = 0.0158;
    b = 0.01084;
    K = 0.01173;
    R = 0.04818;
    L = 0.0629;
else
    load estParameters
    J = estParameters(1).value;
    K = estParameters(2).value;
    L = estParameters(3).value;
    R = estParameters(4).value;
    b = estParameters(5).value;
end
exp_data = ParameterEstimator.TransientExperiment('DC_Motor_ident_m');

% Assign input data.
set(exp_data.InputData(1), 'Data', In(1).values, 'Time', t);

% Assign output data.
set(exp_data.OutputData(1), 'Data', Out(1).values, 'Time', t);
set(exp_data.OutputData(2), 'Data', Out(2).values, 'Time', t);

%% Ident - Set paramteres to estimate
est = ParameterEstimator.Estimation('DC_Motor_ident_m');
est.OptimOptions.Method = 'fmincon';
est.OptimOptions.Algorithm = 'active-set';
est.OptimOptions.TolFun = 1e-4;
est.OptimOptions.TolX = 1e-4;

% Select J
est.Parameters(1).Estimated = true;
set(est.Parameters(1), 'Minimum', 0.0000001, 'Maximum', 1);

% Select K
est.Parameters(2).Estimated = true;
set(est.Parameters(2), 'Minimum', 0.001, 'Maximum', 1);

% Select L
est.Parameters(3).Estimated = true;
set(est.Parameters(3), 'Minimum', 0.0000001, 'Maximum', 0.01);

% Select R
est.Parameters(4).Estimated = true;
set(est.Parameters(4), 'Minimum', 0.05, 'Maximum', 10);

% Select b
est.Parameters(5).Estimated = true;
set(est.Parameters(5), 'Minimum', 0.000001, 'Maximum', 20);

%% Specify the experimental data for the estimation:
est.Experiments = exp_data;

%% Display the estimation iterations.
est.OptimOptions.Display ='iter';

%% Estimate parameters and state.
est.estimate

%% Plot measured data and final simulation responses.
[T,X,Y] = sim('DC_Motor_ident_m', t, [], [t In(1).values]);

figure
subplot(2,1,1)
    plot(t, Out(1).values , T, Y(:,1), '-');
    legend( 'Measured','Simulated');
subplot(2,1,2)
    plot(t, Out(2).values , T, Y(:,2), '-');
    legend( 'Measured','Simulated');

J = est.Parameters(1).value;
K = est.Parameters(2).value;
L = est.Parameters(3).value;
R = est.Parameters(4).value;
b = est.Parameters(5).value;

A = [0      1       0
     0      -b/J    K/J
     0      -K/L    -R/L];

B = [0;...
     0;
     1/L];
C = [1   0  0;
     0   1  0];
 
%% Save paramteres
estParameters = est.Parameters;
save estParameters estParameters

And This is the results:

The paramters after optimization are:

J = 6.9347e-005;
b = 2.2480e-004;
K = 0.0094;
R = 0.2124;
L = 0.0015;

Now the question is why I need to do this in Maple? In my opinion Maple is very  comfortable for analyzing ODE equation. 

Now the goal is to find parameters J, b, K, R, and L in similar way like in Matlab.

@Carl Love Thanks for your comments. V is not function. Is input vector (voltage applied to DC motor) from csv file.

@Preben Alsholm Your interpretation  is  correct: "input signal" =   V(t) . "position of rotor" = theta(t), "speed of rotor"  = diff(theta(t),t). However, keep it mind that V is not a function V(t) is input vector from csv file.

@Preben Alsholm It was a example... Ok. I close this discussion and maybe I will find the solution. Thanks for Your time and Your comments and I apologise for not clear explanation of problem.

@Preben Alsholm  Column 1 - time, column 2 - input signal, column 3 - position of rotor, column 4 - speed of rotor. V(10) = 11.345 this was only for example. Thanks for code.

I found a one solution to my problem. We cane use the following apprche: 

restart;
dsys:={diff(x(t),t)=V,diff(y(t),t)=-x(t),x(0)=1,y(0)=0};
dsn1:=dsolve(dsys,numeric,parameters=[V]);

# Set input value
dsn1(parameters=[10]):
# Get response at time = 1[s]
dsn1(1)

Now using for loop we can find the response of ode system for specific input signal. However, using a loop slows down the program ...

@tomleslie Thank you for your answer, however, I need respons of my system for specific time and input vector.
If I have sinusidal input and length of this signal is 100 and time vector have also length equal to 100 then response of ode also will be of length 100.

@Preben Alsholm From data I need to estimate  J, b, K, R, L.  Maybe I did not clear explain what I want to do. We have ode equation. The input is V - this is voltage applied to dc motor and we have time.  Now for example V(10[s]) = 11.345 volts so I need responso of my ode at time 10[s] and for V = 11.345 volts.  Next  V(11[s]) = 12.345 volts so I need responso of my ode at time 11[s] and for V = 11.345 volts. etc. Then if I have sinusoidal input, my response of ode will be also sinusoidal.My problem is similar to curve fitting however.

For above problem we can use DiffEquation function and use Response Plot or Simulate function but I would like to use  code form my first post.

Below my data: Column 1 - time, column 2 - input signal, column 3 - position of rotor, column 4 - speed of rotor

servo_open_ident_1_A.zip

@Preben Alsholm Yes, each element in the time vector corresponds to an element in the vector V so We can se that V is funtion fo t. For example in csv file I have: time vector, input vector and response of real system. I want to use this data and ode equation for parameter estimation.