Matlab to Simulink optimisation using Ga

Hello,

I want to run simulink model in a function.m for minimisation of the parameter using GA optimisation. so the problem is that the variable output in simulink can be defined in workspace. Then the problem is like this:

The main function is:

clear all;
clc;
pop=1;
Pdim1=100+(600-100)*rand(1,pop);%Pu
Rrl1=0.1+(2.5-0.1)*rand(1,pop);%Rrl
Rdr1=0.3+(0.5-0.03)*rand(1,pop);%Rrl
By1=1.2+(2.2-1.2)*rand(1,pop);%By
p1= 1+(20-1)*rand(1,pop);%p
Js1= 0.1+(5-0.1)*rand(1,pop);%Js
Nepp1= 1+(5-1)*rand(1,pop);%Nepp
vdim1= 3+(32-3)*rand(1,pop);%Vitesse
Udim1= (150)*rand(1,pop)%Tension
pmin=[100 0.1 0.03 1.2 1 0.1 1 3 0];
pmax=[600 2.5 0.5 2.2 20 5 5 32 150];
params=[Pdim1 Rrl1 Rdr1 By1 p1 Js1 Nepp1 vdim1 Udim1];
options=gaoptimset(display,'iter');
options=gaoptimset(Generations,'500')
[y,x]=ga(@objectif,9,[],[],[],[],pmin,pmax,[],[]);
[popt,Fx]=objectif(params);

and the objectif function is:

function [popt,Fx]=objectif(params);
%saisi de tout les constante de la machine:
Kc=1; %Coefficient de carter
Kr=0.5; % Coefficient de remplissage des encoches de la génératrice
a_aimant=pi/3; %angle de l aimant
Br=1.1;
Mr=1.05; %permeablite relative de l'aimant 
Kb=0.95; %Coefficient de bobinage
Mo=4*pi*10^(-7); %Permeabilite du vide
ro_cuivre=8960 % la masse volumique du cuivre
ro_fer=7860 % la masse volumique du matériau ferromagnétique de la culasse
ro_aimant=3450 % la masse volumique daimant
Kp=0.833;
K_r=0.35;
lmg=4;

%------------------------------ Condition et limite des variables -----------------------------%
% Pu=[100,600]; Rapport rayon d'alésage/longueur=[0.1,2.5]; Induction dans la culasse=[1.2,2.2];
% p=[1,20]; densite de courant dans les encoches=[0.1,5]; Nombre d'encoche par pole par phase=[1,5];
% vitesse de rotation=[3,32]; Tension =[0,150];

%%

%------------------------------ Intialization -----------------------------%
%%

%
Pdim1=params(1);%Pu
Rrl1=params(2);%Rrl
Rdr1=params(3);%Rrl
By1=params(4);%By
p1= params(5);%p
Js1= params(6);%Js
Nepp1= params(7);%Nepp
vdim1= params(8);%Vitesse
Udim1= params(9)%Tension
   
Pu=Pdim1;
Rrl=Rrl1;
Rdr=Rdr1;
By=By1;
p= p1;
Js= Js1;
Nepp= Nepp1;
vdim= vdim1;
Udim= Udim1;


rs=0.025; %on va prendre comme constante la valeur de rs (rayon d allesage)
lr=rs/Rrl;
ds=rs*Rdr;
Nenc=6*p*round(Nepp);
wT=(4*pi*rs)/(3*(round(Nepp)));
ws=wT;
g=0.001+(0.003*sqrt(rs*lr));
h2=ws/8;
h3=0.02*rs;
h1=8*ds*Kr/7;
b1=ws;
b2=ws/2;
b3=3*ws/4;
Ba=Br*(lmg/(Mr+lmg));
B_1a=(4/pi)*Ba*sin(a_aimant);
lm=Kc*g*(Mr/((Br/Ba)-1));
dr=(rs*Ba*a_aimant)/(p*By);
dy=dr;
wm=(2*a_aimant*rs)/p;
S_enc=0.5*(b1+b3)*h1*Kr;

K_z1=sin(pi/6)/(Nepp*sin(pi/(6*Nepp)));
K_1b=K_z1;
K_1s=Js*ws*ds*K_1b*Kr/(ws+wT);
Tem=pi*rs^2*lr*B_1a*Js*ds*Kr*K_1b;

%%
%Calcule de parametre electromangetique et mecanique de la generatrice
%%
Fi_aimant=Ba*lr*rs*a_aimant; %Flux maximale produite par l aimant
Fi_y=Fi_aimant/2; %Flux canaliser par la culasse
Fi_s1=2*K_1b*round(Nepp)*B_1a*rs*lr; %Flux totale
Fi_eff1=Fi_s1;
Lm1=(4*Mo*lr*rs*(K_1b)^2*(round(Nepp))^2)/(pi*(Kc*g+(lm/Mr))); % Inductance magnetisante
l_enc=(2*h1/(3*(b1+b3)))+(2*h2/(b2+b3))+(h3/b2);
Lf1=2*Mo*lr*p*Nepp*l_enc; %Inductance de fuite
M1=-Lm1/2;
L1=Lm1-M1+Lf1;%inductance totale
I1=Fi_s1/L1;
l_tete=pi*(rs+(ds/2))/p;
sigma=59.6*10^6; %Conductivite thermik du cuivre
R1=2*p*(lr+l_tete)*Nepp/(sigma*S_enc);
% A3=1;
% be2=(2*Udim*R1*I1)/(((Fi_s1*p*vdim)^2-(R1^2+(L1*p*vdim)^2))*I1^2);
% B3=Udim^2/(((Fi_s1*p*vdim)^2-(R1^2+(L1*p*vdim)^2))*I1^2);
% if ((be2^2-B3)>0)
%     Nce=sqrt(be2^2-B3)-be2;
% else
    Nce=120;
% end


Rs=R1*Nce^2;
Ls=L1*Nce^2;
Fi_eff=Fi_s1*Nce
%------------------------- Calcule de la masse en general ------------------------%

%%
Rint=rs-g-lm-dr;
rrotor=rs-g-lm;
Ra=rs-g;%ws Largeur de dent);
Mcr=pi*lr*(rrotor^2-Rint^2)*ro_fer;
Rint=rs-g-lm-dr;%rayon interieur de l aimant
rrotor = rs-g-lm;
V_aimant=pi.*lr.*(Kp.*(Ra.^2-rrotor.^2));
M_aimant=ro_aimant.*V_aimant;
M_rotor=Mcr+M_aimant;
Rcs=rs+ds;
Rext=rs+ds+dy;
Vcs=2*lr*dy*(rs+ds+dy/2);
Mcs=ro_fer*Vcs;
Teta_b=(2*pi./Nenc)-(b2/rs);
V_Pdent=Nenc.*rs.*Teta_b.*(2.*h2+h3).*(lr/2); %volume des pieds de dents
V_hdent=pi.*lr.*((rs.*(ds-(h2+h3)))+((ds.^2-(h2+h3).^2))/2);% volume des dents hors isthme 
M_dents=ro_fer*(V_Pdent+V_hdent);
M_stator=Mcs+M_dents;
Vcuivre_enc=pi*lr*ds*Kr*(ds+(rs/2));
ltetes_bob=pi^2*(rs+0.5*ds)/(2*p);
Vtete_bob=pi*ltetes_bob*ds*ws*Nenc*Kr;
M_cuivre=(Vcuivre_enc+Vtete_bob)*ro_cuivre;
M_gen=M_rotor+M_stator+M_cuivre;

lambda=0.75; %La valeur de TSR tip speed rotor
r=0.5; %rayon de l'éolienne
S=pi*r*1.5; %surface ballayer par le vent 
ro=1.292; %masse volumique de l'air
Jturb=0.01; %*M_gen*rs^2; %Couple de la turbine
Jgen=(0.5*M_rotor*(rs+dy)^2)*1e3;%Inertie de la génératrice
Jtur=Jturb+Jgen;% Tgen=5 %couple de la génératrice 
G=3; %multiplicateur de vitesse
Fmec=0.06 ;%frottement 
Rch=0;
Rtotal=Rs+Rch;
Lch=0;
assignin('base','Rs',Rs);
assignin('base','Ls',Ls);
assignin('base','p',p);
assignin('base','Fi_eff',Fi_eff);
assignin('base','Jturb',Jturb);
assignin('base','Lch',Lch);
assignin('base','Rch',Rch);
assignin('base','G',G);
assignin('base','Fmec',Fmec);
open('MSAP2');
sim('MSAP2','SrcWorkspace','Current');
open('RNX');
sim('RNX','SrcWorkspace','Current');
global Pdc
get Pdc;
pmaxi=norm(Pdc);
Fx=-pmaxi;
popt=min(Fx);

Answer 1

[Sulaymon Eshkabilov]

Hi

If I've understood your question correctly is that you are trying to run/simulate your pre-existing Simulink model from your above shown m-file, after initialinzing some input variable data for your simulink model, right? And then to run ga based optimization to get optimal values to be fed back your Simulink model, correct?

If so, then you'd need to employ the commands: sim() to run your simulink model within your m-file and get_param() to obtain the simulation parameters of your simulink model.

That is all.

Good luck

[RANJORO Victor Nantenaina]

Yes and the problem is that the variable under the function still in and it can be shown in workspace. I want to optimize the output 'Pdc' where the variable 'params' using ga for input. these are the program main and the objectif function