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);%RrlBy1=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);%PuRrl1=params(2);%RrlRdr1=params(3);%RrlBy1=params(4);%Byp1= params(5);%pJs1= params(6);%JsNepp1= params(7);%Neppvdim1= params(8);%VitesseUdim1= 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 Pdcget Pdc;pmaxi=norm(Pdc);Fx=-pmaxi;popt=min(Fx);
Best Answer