function R = resistance_fissure_2D(Lf, ef)
%% Dimensions en mm
L = 100; % longueur de la barre [mm]
H = 50; % largeur dans le plan 2D [mm]
p = 30; % profondeur [mm]
%% Donnees electriques
sigma = 5.8e7; % conductivite [S/m]
Vg = 1; % potentiel gauche [V]
Vd = 0; % potentiel droite [V]
%% Fissure
% Lf : hauteur de fissure [m]
% ef : largeur de fissure [m]
xc = L/2;
x1f = xc - ef/2;
x2f = xc + ef/2;
y1f = H - Lf;
y2f = H;
%% Ouverture FEMM
openfemm(1);
newdocument(3); % 3 = current flow problem
%% Definition du probleme
ci_probdef('millimeters','planar',0,1e-8,p,30);
%% Materiaux
ci_addmaterial('Metal', sigma, sigma, 1, 1, 0, 0);
%% Conducteurs
ci_addconductorprop('Gauche', Vg, 0, 1);
ci_addconductorprop('Droite', Vd, 0, 1);
%% Geometrie de la barre
ci_addnode(0,0);
ci_addnode(L,0);
ci_addnode(L,H);
ci_addnode(0,H);
ci_addsegment(0,0,L,0);
ci_addsegment(L,0,L,H);
ci_addsegment(0,H,0,0);
%% Geometrie de la fissure
ci_addnode(x1f,y1f);
ci_addnode(x2f,y1f);
ci_addnode(x2f,y2f);
ci_addnode(x1f,y2f);
ci_addsegment(x1f,y1f,x2f,y1f);
ci_addsegment(x2f,y1f,x2f,y2f);
ci_addsegment(x1f,y2f,x1f,y1f);
ci_addsegment(x1f,y2f,0,H);
ci_addsegment(x2f,y2f,L,H);
%% Affectation des conducteurs aux bords gauche et droit
ci_clearselected;
ci_selectsegment(0,H/2);
ci_setsegmentprop('<None>',0,1,0,0,'Gauche');
ci_clearselected;
ci_selectsegment(L,H/2);
ci_setsegmentprop('<None>',0,1,0,0,'Droite');
%% Bloc conducteur
ci_addblocklabel(L/4,H/2);
ci_selectlabel(L/4,H/2);
ci_setblockprop('Metal',1,0,1);
ci_clearselected;
%% La fissure est vide : pas de materiau dedans
% On ne met pas de block label dans la fissure.
%% Sauvegarde et calcul
nomfichier = fullfile(pwd,'barre_fissuree.FEC');
ci_saveas(nomfichier);
ci_analyze(1);
ci_loadsolution;
%% Recuperation du courant
Val_gauche = co_getconductorproperties('Gauche');
Val_droite = co_getconductorproperties('Droite');
I = Val_gauche(2);
R = abs(Val_gauche(1) - Val_droite(1))/I;
%% Fermeture
co_close;
ci_close;
end
clear; clc; close all;
%% Fissure rÈelle inconnue
Hf_reel = 18; % c'est la valeur reelle de la hauteur du defaut (supposee connue)
Lf_fixe = 2; % C'est la valeur reelle de la largeur du defaut qui est a priori connue
Rmes = resistance_fissure_2D(Hf_reel, Lf_fixe); % on calcul la resistance aves les
valeurs reelles
fprintf('RÈsistance mesurÈe simulÈe = %.6e Ohm\n', Rmes);
%% Bornes sur Hf
Hf_min = 2;
Hf_max = 30;
%% Recherche globale par grille 1D
Hf_vec = linspace(Hf_min, Hf_max, 50);
Jvec = zeros(length(Hf_vec),1);
Rvec = zeros(length(Hf_vec),1);
bestJ = inf;
bestHf = NaN;
for j = 1:length(Hf_vec)
Hf = Hf_vec(j);
R = resistance_fissure_2D(Hf, Lf_fixe);
J = ((R - Rmes)/Rmes)^2;
Rvec(j) = R;
Jvec(j) = J;
if J < bestJ
bestJ = J;
bestHf = Hf;
end
end
fprintf('\nMeilleur point grille :\n');
fprintf('Hf = %.3f mm\n', bestHf);
fprintf('Lf fixÈ = %.3f mm\n', Lf_fixe);
fprintf('J = %.6e\n', bestJ);
%% Raffinement local avec fminsearch
Jfun = @(Hf) objectif_fissure_1D(Hf, Rmes, Lf_fixe, Hf_min, Hf_max);
options = optimset('Display','iter', ...
'TolX',1e-9, ...
'TolFun',1e-14, ...
'MaxIter',100);
Hf_opt = fminsearch(Jfun, bestHf, options);
Ropt = resistance_fissure_2D(Hf_opt, Lf_fixe);
fprintf('\nRÈsultats identification 1 paramËtre :\n');
fprintf('Hf rÈel = %.3f mm\n', Hf_reel);
fprintf('Hf identifiÈ = %.3f mm\n', Hf_opt);
fprintf('Lf fixÈ = %.3f mm\n', Lf_fixe);
fprintf('R mesurÈe = %.6e Ohm\n', Rmes);
fprintf('R identifiÈe = %.6e Ohm\n', Ropt);
fprintf('Erreur rel. R = %.6f %%\n', abs(Ropt-Rmes)/Rmes*100);
%% Affichage J(Hf)
figure;
plot(Hf_vec*1e3, log10(Jvec + 1e-16), 'b-', 'LineWidth', 2);
hold on;
plot(Hf_reel*1e3, log10(1e-16), 'wo', 'MarkerSize', 10, 'LineWidth', 2);
plot(Hf_opt*1e3, log10(Jfun(Hf_opt) + 1e-16), 'rx', 'MarkerSize', 10, 'LineWidth', 2);
grid on;
xlabel('Longueur fissure L_f [mm]');
ylabel('log_{10}(J)');
title('Identification de L_f avec e_f fixÈ');
legend('J(L_f)', 'Fissure rÈelle', 'Fissure identifiÈe');
%% Affichage R(Hf)
figure;
plot(Hf_vec, Rvec, 'k-', 'LineWidth', 2);
hold on;
plot(Hf_reel, Rmes, 'wo', 'MarkerSize', 10, 'LineWidth', 2);
plot(Hf_opt, Ropt, 'rx', 'MarkerSize', 10, 'LineWidth', 2);
grid on;
xlabel('Longueur fissure L_f [mm]');
ylabel('RÈsistance R [Ohm]');
title('RÈsistance Èlectrique en fonction de L_f');
legend('R(L_f)', 'Fissure rÈelle', 'Fissure identifiÈe');
%% Fonction objectif 1D
function val = objectif_fissure_1D(Hf, Rmes, Lf_fixe, Hf_min, Hf_max)
if Hf < Hf_min || Hf > Hf_max
val = 1e6;
return;
end
R = resistance_fissure_2D(Hf, Lf_fixe);
val = ((R - Rmes)/Rmes)^2;
end
clear; clc; close all;
%% Fissure rÈelle
Hf_reel = 18; % C'est la valeur reelle de la hauteur du defaut (supposee connue)
Lf_fixe = 2; % C'est la valeur reelle de la largeur du defaut qui est a priori connue
Rmes = resistance_fissure_2D(Hf_reel, Lf_fixe); % on calcul la resistance aves les
valeurs reelles
fprintf('RÈsistance mesurÈe simulÈe = %.6e Ohm\n', Rmes);
%% Bornes de recherche pour Hf
Hf_min = 2;
Hf_max = 30;
lb = Hf_min;
ub = Hf_max;
%% Fonction objectif
J = @(Hf) objectif_PSO_1param(Hf, Lf_fixe, Rmes);
%% Options PSO
options = optimoptions('particleswarm', ...
'SwarmSize', 20, ...
'MaxIterations', 30, ...
'Display', 'iter', ...
'FunctionTolerance', 1e-12);
%% Lancement PSO
nvars = 1;
[Hf_opt, Jopt] = particleswarm(J, nvars, lb, ub, options);
Ropt = resistance_fissure_2D(Hf_opt, Lf_fixe);
%% RÈsultats
fprintf('\nRÈsultats identification par PSO, 1 paramËtre :\n');
fprintf('Hf rÈel = %.3f mm\n', Hf_reel);
fprintf('Hf identifiÈ = %.3f mm\n', Hf_opt);
fprintf('Lf fixÈ = %.3f mm\n', Lf_fixe);
fprintf('R mesurÈe = %.6e Ohm\n', Rmes);
fprintf('R identifiÈe = %.6e Ohm\n', Ropt);
fprintf('Erreur rel. R = %.6f %%\n', abs(Ropt-Rmes)/Rmes*100);
fprintf('J optimal = %.6e\n', Jopt);
%% TracÈ de la fonction objectif
Hf_vec = linspace(Hf_min, Hf_max, 60);
J_vec = zeros(size(Hf_vec));
for i = 1:length(Hf_vec)
J_vec(i) = objectif_PSO_1param(Hf_vec(i), Lf_fixe, Rmes);
end
figure;
plot(Hf_vec, J_vec, 'LineWidth', 2);
hold on;
plot(Hf_reel, 0, 'ro', 'MarkerSize', 8, 'LineWidth', 2);
plot(Hf_opt, Jopt, 'kx', 'MarkerSize', 10, 'LineWidth', 2);
grid on;
xlabel('Longueur de fissure L_f [mm]');
ylabel('Fonction objectif J');
title('Identification de L_f par PSO avec e_f fixÈ');
legend('J(L_f)', 'Valeur rÈelle', 'Valeur identifiÈe');
%% Fonction objectif locale
function J = objectif_PSO_1param(Hf, Lf_fixe, Rmes)
R = resistance_fissure_2D(Hf, Lf_fixe);
J = ((R - Rmes)/Rmes)^2;
end
clear; clc; close all;
%% Fissure rÈelle
Hf_reel = 18; % c'est la valeur reelle de la hauteur du defaut (supposee connue)
Lf_fixe = 2; % C'est la valeur reelle de la largeur du defaut qui est a priori connue
Rmes = resistance_fissure_2D(Hf_reel, Lf_fixe); % on calcul la resistance aves les
valeurs reelles
fprintf('RÈsistance mesurÈe simulÈe = %.6e Ohm\n', Rmes);
%% Bornes de Hf
Hf_min = 2;
Hf_max = 30;
%% ParamËtres EGO
n_init = 6;
n_iter = 20;
n_cand = 500;
%% Points initiaux
X = Hf_min + rand(n_init,1)*(Hf_max-Hf_min);
Y = zeros(n_init,1);
for i = 1:n_init
Y(i) = objectif_1param(X(i), Lf_fixe, Rmes);
end
%% Boucle EGO
for iter = 1:n_iter
[Ybest, idx] = min(Y);
xbest = X(idx);
fprintf('ItÈration %d : Jbest = %.6e, Hf = %.3f mm\n', ...
iter, Ybest, xbest);
Xcand = Hf_min + rand(n_cand,1)*(Hf_max-Hf_min);
[mu, s] = modele_RBF_1D(X, Y, Xcand);
EI = expected_improvement(Ybest, mu, s);
[~, ind] = max(EI);
xnew = Xcand(ind);
ynew = objectif_1param(xnew, Lf_fixe, Rmes);
X = [X; xnew];
Y = [Y; ynew];
end
%% RÈsultat final
[Yopt, idx] = min(Y);
Hf_opt = X(idx);
Ropt = resistance_fissure_2D(Hf_opt, Lf_fixe);
fprintf('\nRÈsultats identification EGO 1 paramËtre :\n');
fprintf('Hf rÈel = %.3f mm\n', Hf_reel);
fprintf('Hf identifiÈ = %.3f mm\n', Hf_opt);
fprintf('Lf fixÈ = %.3f mm\n', Lf_fixe);
fprintf('R mesurÈe = %.6e Ohm\n', Rmes);
fprintf('R identifiÈe = %.6e Ohm\n', Ropt);
fprintf('Erreur rel. R = %.6f %%\n', abs(Ropt-Rmes)/Rmes*100);
fprintf('J optimal = %.6e\n', Yopt);
%% TracÈ J(Hf)
Hf_plot = linspace(Hf_min, Hf_max, 60)';
J_plot = zeros(size(Hf_plot));
for i = 1:length(Hf_plot)
J_plot(i) = objectif_1param(Hf_plot(i), Lf_fixe, Rmes);
end
figure;
plot(Hf_plot, J_plot, 'LineWidth', 2);
hold on;
plot(Hf_reel, 0, 'ro', 'MarkerSize', 8, 'LineWidth', 2);
plot(Hf_opt, Yopt, 'kx', 'MarkerSize', 10, 'LineWidth', 2);
grid on;
xlabel('Longueur fissure L_f [mm]');
ylabel('Fonction objectif J');
title('Identification de L_f avec e_f fixÈ');
legend('J(L_f)','Valeur rÈelle','Valeur identifiÈe');
%% Fonctions locales
function J = objectif_1param(Hf, Lf_fixe, Rmes)
R = resistance_fissure_2D(Hf, Lf_fixe);
J = ((R - Rmes)/Rmes)^2;
end
function [mu, s] = modele_RBF_1D(X, Y, Xcand)
n = length(X);
m = length(Xcand);
longueur = 5e-3;
D = zeros(n,n);
for i = 1:n
for j = 1:n
r = abs(X(i)-X(j));
D(i,j) = exp(-(r/longueur)^2);
end
end
lambda = 1e-10;
coef = (D + lambda*eye(n)) \ Y;
mu = zeros(m,1);
s = zeros(m,1);
for k = 1:m
d = zeros(n,1);
for i = 1:n
r = abs(Xcand(k)-X(i));
d(i) = exp(-(r/longueur)^2);
end
mu(k) = d' * coef;
dist_min = min(abs(X - Xcand(k)));
s(k) = dist_min;
end
s = s / max(s + eps);
s = s + 1e-12;
end
function EI = expected_improvement(Ybest, mu, s)
z = (Ybest - mu)./s;
Phi = 0.5*(1 + erf(z/sqrt(2)));
phi = exp(-0.5*z.^2)/sqrt(2*pi);
EI = (Ybest - mu).*Phi + s.*phi;
EI(EI < 0) = 0;
end