Laboratorul 1 – MATLAB (1)

Recapitulare: vectori, matrici, operatori logici și reprezentări grafice

Descarcă codul MATLAB (.m)

Lucrarea recapitulează lucrul cu vectori și matrici în MATLAB, operațiile element-wise, operatorii relaționali și logici, selecția elementelor cu indecși logici sau cu funcția find, ciclurile for și reprezentările grafice 2D/3D.

Cum rulezi codul din laborator3 min

Film în pregătire

Deschizi fișierul cu cod în MATLAB, rulezi o singură secțiune cu Run Section și urmărești variabilele în Workspace.

Deschide filmul pe pagina lui

Competențe urmărite

definirea vectorilor și generarea automată a valorilor;
calculul tabelar folosind operații vectoriale și matriciale;
folosirea operatorilor relaționali și a operatorilor logici &, |, ~;
selectarea datelor pe baza unor condiții logice;
reorganizarea unui vector cu for și cu indexare logică;
reprezentarea grafică 2D și 3D a rezultatelor;
interpretarea rezultatului numeric în contextul unei aplicații inginerești.

1. Noțiuni teoretice MATLAB

MATLAB lucrează natural cu vectori și matrici. Într-un laborator de modelare și simulare, această proprietate este esențială: o singură expresie poate calcula simultan valori pentru mai multe momente de timp, mai mulți parametri sau mai multe scenarii de funcționare.

1.1 Vectori și matrici

Un vector linie se poate nota sub forma

x = [\,x_1 \;\; x_2 \;\; \dots \;\; x_n\,]

iar un vector coloană se obține prin transpunere. În MATLAB, transpunerea se realizează cu apostrof:

x = [1 2 3 4 5];
y = x';

O matrice conține elemente dispuse pe linii și coloane. Elementul de pe linia i și coloana j se notează uzual cu a_{ij}. În MATLAB, accesarea acestuia se face cu A(i,j).

1.2 Operații element-wise

În MATLAB trebuie făcută diferența între operațiile matriciale și operațiile aplicate element cu element (element-wise). Pentru vectori și semnale eșantionate, cele mai frecvente operații sunt element-wise.

Operator	Semnificație	Exemplu
`.*`	înmulțire element cu element	`y = a.*b`
`./`	împărțire element cu element	`y = a./b`
`.^`	putere element cu element	`y = x.^2`
`*`	produs matricial	`C = A*B`
`^`	putere matricială	`A^2`

Atenție

Pentru un vector x, expresia x^2 este în general eronată, deoarece MATLAB încearcă o operație matricială. Pentru ridicarea la pătrat a fiecărui element se folosește x.^2.

Operații element cu element: x.^2 și x^23 min

Film în pregătire

Ce eroare dă x^2 pe un vector, cum o citești și de ce x.^2 este varianta corectă. La fel pentru .* și ./.

Deschide filmul pe pagina lui

1.3 Operatori relaționali, logici și funcția find

Operatorii relaționali produc vectori logici. Operatorii logici permit combinarea mai multor condiții. În problemele inginerești, aceștia sunt utili pentru identificarea intervalelor sigure de funcționare, a depășirilor de prag sau a valorilor invalide dintr-un semnal.

Operator	Semnificație	Exemplu
`<`	mai mic	`x < 0`
`<=`	mai mic sau egal	`T <= 60`
`>`	mai mare	`v > 8`
`>=`	mai mare sau egal	`h >= 12`
`==`	egal	`k == 3`
`~=`	diferit	`x ~= 0`
`&`	ȘI logic	`h >= 12 & abs(v) <= 8`
`\|`	SAU logic	`B < 30 \| Temp > 70`
`~`	negație	`~masca`

Funcția find returnează pozițiile elementelor nenule sau ale elementelor pentru care o condiție logică este adevărată:

idx = find(x > 0);
valori_pozitive = x(idx);

Problema 1.1 – Energia cinetică a unui vehicul

Energia cinetică a unui vehicul de masă m, care se deplasează cu viteza v, este:

E_k = \tfrac{1}{2}\, m\, v^2

Se consideră masele m = 800, 1000, 1200, 1400, 1600 kg și vitezele v = 20, 40, 60, 80, 100 km/h. Să se creeze un tabel în care liniile corespund maselor, iar coloanele corespund vitezelor. Valorile energiei se vor exprima în kJ.

Observație

Viteza trebuie convertită din km/h în m/s înainte de aplicarea formulei energiei cinetice, deoarece unitățile din Sistemul Internațional sunt kg, m și s.

%% Problema 1.1 - Energie cinetică pentru un vehicul electric
clear; clc; close all

m = 800:200:1600;        % mase [kg]
v_kmh = 20:20:100;       % viteze [km/h]
v_ms = v_kmh/3.6;        % conversie în [m/s]

Ek = 0.5*(m')*(v_ms.^2); % energie [J]
Ek_kJ = Ek/1000;         % energie [kJ]

T_Ek = array2table(Ek_kJ, ...
    'VariableNames', compose('v_%d_kmh', v_kmh), ...
    'RowNames', compose('m_%d_kg', m));

disp(T_Ek)

figure
plot(v_kmh, Ek_kJ', 'LineWidth', 1.5)
grid on
xlabel('Viteza [km/h]')
ylabel('Energia cinetică [kJ]')
title('Energia cinetică pentru mase diferite')
legend(compose('m = %d kg', m), 'Location', 'northwest')

Problema 1.1: un tabel dintr-un singur produs4 min

Film în pregătire

Cum produsul dintre un vector coloană și un vector linie dă toată matricea energiilor, apoi tabelul cu array2table și graficul.

Deschide filmul pe pagina lui

Rezultatul rulării

Tabelul 1.1. Energia cinetică [kJ] pentru combinațiile masă–viteză analizate.
Masa [kg]	20 km/h	40 km/h	60 km/h	80 km/h	100 km/h
800	12.35	49.38	111.11	197.53	308.64
1000	15.43	61.73	138.89	246.91	385.80
1200	18.52	74.07	166.67	296.30	462.96
1400	21.60	86.42	194.44	345.68	540.12
1600	24.69	98.77	222.22	395.06	617.28

Interpretare. Pentru aceeași masă, energia cinetică crește cu pătratul vitezei. De exemplu, dublarea vitezei de la 20 km/h la 40 km/h conduce la o energie cinetică de patru ori mai mare.

Energia cinetică în funcție de viteză, câte o curbă pentru fiecare masă

În expresia (m')*(v_ms.^2), primul vector este coloană, iar al doilea este linie. Rezultatul este o matrice cu 5 linii și 5 coloane. Operatorul .^ este obligatoriu deoarece ridicarea la pătrat se aplică fiecărui element al vectorului vitezelor.

Problema 1.2 – Fereastra sigură de urcare pentru o dronă

O dronă execută o urcare verticală comandată. Pentru descrierea nominală a altitudinii se consideră o viteză medie de urcare v_c constantă. Măsurarea vitezei verticale prezintă însă o variație sinusoidală, care modelează simplificat efectul unei perturbații verticale (de exemplu, rafale de vânt).

h(t) = h_0 + v_c\, t

v_z(t) = v_c + A \sin(\omega t) \qquad \text{(viteza verticală reală)}

Se consideră h_0 = 2 m, v_c = 3 m/s, A = 2 m/s și \omega = 0,8 rad/s, pentru 0 ≤ t ≤ 10 s. Determinați intervalul de timp în care drona se află simultan între 12 m și 25 m altitudine, iar viteza verticală măsurată nu depășește 3,5 m/s. Folosiți operatori logici și relaționali și determinați intervalul atât cu funcția find, cât și prin indexare logică.

%% Problema 1.2 - Fereastra sigură de urcare pentru o dronă
clear; clc; close all

h0 = 2;          % altitudine inițială [m]
vc = 3;          % viteză verticală medie [m/s]
A = 2;           % amplitudinea perturbației [m/s]
omega = 0.8;     % pulsație [rad/s]

t = linspace(0,10,501);
h = h0 + vc*t;                         % model nominal al altitudinii
vz = vc + A*sin(omega*t);              % viteză verticală măsurată

% Condițiile ferestrei sigure de urcare
conditie = (h >= 12) & (h <= 25) & (vz <= 3.5);
idx = find(conditie);

t1 = t(idx(1));
t2 = t(idx(end));
fprintf('Metoda 1: interval aproximativ [%.2f, %.2f] s\n', t1, t2)

% Metoda 2: indexare logică
t_selectat = t(conditie);
fprintf('Metoda 2: interval aproximativ [%.2f, %.2f] s\n', ...
    t_selectat(1), t_selectat(end))

figure
plot(t,h,'LineWidth',1.5); hold on
plot(t(conditie),h(conditie),'LineWidth',3)
yline(12,'--'); yline(25,'--')
grid on
xlabel('t [s]'); ylabel('h [m]')
title('Fereastra sigură - condiția de altitudine')

figure
plot(t,vz,'LineWidth',1.5); hold on
plot(t(conditie),vz(conditie),'LineWidth',3)
yline(3.5,'--')
grid on
xlabel('t [s]'); ylabel('v_z [m/s]')
title('Fereastra sigură - condiția de viteză')

Problema 1.2: condiții logice, find și indexare logică4 min

Film în pregătire

Cum arată vectorul logic conditie în Workspace, cum se obține intervalul cu find și cu indexare logică și cum apare porțiunea selectată pe grafic.

Deschide filmul pe pagina lui

Altitudinea dronei în timp, cu porțiunea din fereastra sigură îngroșată și pragurile de 12 m și 25 m

Viteza verticală a dronei în timp, cu porțiunea din fereastra sigură îngroșată și pragul de 3,5 m/s

Condiția compusă utilizează operatorul logic AND (&): toate cerințele trebuie îndeplinite simultan. În acest exemplu, altitudinea este descrisă printr-un model nominal simplu, h(t) = h_0 + v_c t, iar variația sinusoidală este introdusă în semnalul vitezei verticale măsurate pentru a modela apariția unor perturbații (rafale de vânt).

Modelul este intenționat simplificat, deoarece obiectivul problemei este utilizarea operatorilor relaționali, a funcției find și a indexării logice, nu modelarea dinamică detaliată a dronei.

PRO 1.3 – Alarmă pentru baterie și temperatură

Pentru o dronă aflată în test, nivelul de încărcare a bateriei B(t) scade aproximativ liniar, iar temperatura motorului T(t) crește către o valoare de regim staționar. Se consideră:

B(t) = 100 - 4.5\, t

T(t) = T_a + \Delta T \left(1 - e^{-t/\tau}\right)

unde T_a = 25 °C, \Delta T = 60 °C și \tau = 6 min, pentru 0 ≤ t ≤ 20 min. Definiți separat alarma de baterie scăzută și alarma de temperatură crescută, apoi construiți alarma generală, care se activează dacă bateria scade sub 30% SAU temperatura motorului depășește 70 °C. Determinați momentul primei activări și reprezentați grafic cele două cauze ale alarmei și starea alarmei generale.

%% PRO 1.3 - Alarmă pentru baterie și temperatură
clear; clc; close all
t = linspace(0,20,401);       % timp [min]
B = 100 - 4.5*t;              % baterie [%]
Ta = 25;                      % temperatura ambiantă [°C]
DeltaT = 60;                  % creșterea de temperatură la regim stationar [°C]
tau = 6;                      % constanta de timp termică [min]
Temp = Ta + DeltaT*(1-exp(-t/tau));

% Măști logice separate
alarma_B = B < 30;
alarma_T = Temp > 70;

% Alarma generală: este suficientă oricare dintre cele două cauze (SAU)
alarma = alarma_B | alarma_T;

idx = find(alarma);
fprintf('Prima alarmă generală apare la t = %.2f min\n', t(idx(1)))

idxB = find(alarma_B,1,'first');
idxT = find(alarma_T,1,'first');
fprintf('Alarmă temperatură: t = %.2f min\n', t(idxT))
fprintf('Alarmă baterie: t = %.2f min\n', t(idxB))

figure
tiledlayout(3,1)

nexttile
plot(t,B,'LineWidth',1.5); hold on
plot(t(alarma_B),B(alarma_B),'LineWidth',3)
yline(30,'--')
grid on
ylabel('Baterie [%]')
title('Alarmă baterie: B < 30%')

nexttile
plot(t,Temp,'LineWidth',1.5); hold on
plot(t(alarma_T),Temp(alarma_T),'LineWidth',3)
yline(70,'--')
grid on
ylabel('Temperatură [°C]')
title('Alarmă temperatură: T > 70 °C')

nexttile
stairs(t,alarma,'LineWidth',2)
ylim([-0.1 1.1])
yticks([0 1])
yticklabels({'OFF','ON'})
grid on
xlabel('t [min]')
title('Starea alarmei generale: alarma_B | alarma_T')

Trei grafice: nivelul bateriei cu pragul de 30%, temperatura cu pragul de 70 °C și starea alarmei generale

Interpretare. Alarma de temperatură se activează prima, la aproximativ 8,32 min. Alarma de baterie se activează ulterior, la aproximativ 15,56 min. Alarma generală rămâne activă din momentul în care prima dintre cele două condiții devine adevărată.

Problema 1.4 – Curățarea unui semnal experimental

În urma unei măsurători s-a obținut un vector x. Valorile foarte mici, cu modulul mai mic decât 0.05, sunt considerate zgomot neglijabil. Dorim să păstrăm la începutul vectorului valorile utile, în ordinea inițială, iar la sfârșit să mutăm un număr corespunzător de zerouri.

De exemplu: x = [0.82, -0.03, 1.15, 0.04, -2.40, 0.00, 0.51, -0.02, 1.90]

Condiția de prag: |x_i| < 0.05

%% Problema 1.4 - Curățarea unui semnal experimental
clear; clc; close all

x = [0.82, -0.03, 1.15, 0.04, -2.40, 0.00, 0.51, -0.02, 1.90];
prag = 0.05;

% Varianta 1: cu for
utile = [];
eliminate = [];

for k = 1:length(x)
    if abs(x(k)) >= prag
        utile = [utile, x(k)];
    else
        eliminate = [eliminate, x(k)];
    end
end

x_for = [utile, zeros(size(eliminate))]

% Varianta 2: cu indexare logică
masca_utile = abs(x) >= prag;
x_logic = [x(masca_utile), zeros(1,sum(~masca_utile))]

% Verificare
isequal(x_for, x_logic)

Problema 1.4: ciclul for, pas cu pas4 min

Film în pregătire

Rularea ciclului for cu punct de oprire (breakpoint), ca să vezi cum cresc vectorii utile și eliminate, apoi aceeași operație într-o singură linie.

Deschide filmul pe pagina lui

Rezultatul rulării

x_for =
    0.8200    1.1500   -2.4000    0.5100    1.9000         0         0         0         0

x_logic =
    0.8200    1.1500   -2.4000    0.5100    1.9000         0         0         0         0

ans =
  logical
   1

Interpretare. Cele două metode produc același vector. Valoarea logică 1 returnată de isequal confirmă egalitatea rezultatelor.

În practică, varianta cu indexare logică este mai compactă și mai eficientă. Varianta cu for rămâne utilă didactic, deoarece arată explicit cum sunt parcurse elementele vectorului.

Problema 1.5 – Traiectoria 3D a unei drone în prezența vântului

O dronă execută o manevră elicoidală ascendentă în jurul unei axe verticale. În absența vântului, proiecția traiectoriei în planul orizontal este un cerc de rază R, iar altitudinea crește cu viteza verticală constantă v_z. O componentă constantă a vântului pe direcția x (notată cu v_{w,x} mai jos) produce o derivă laterală a întregii traiectorii.

x(t) = R \cos(\omega t) + v_{w,x}\, t

y(t) = R \sin(\omega t)

z(t) = z_0 + v_z\, t

Se consideră R = 12 m, z_0 = 5 m, v_z = 1.5 m/s și \omega = 0.5 rad/s, pentru 0 ≤ t ≤ 25 s. Comparați trei cazuri: fără vânt, vânt spre Vest de 2 m/s și vânt spre Est de 2 m/s. Reprezentați pe același grafic 3D cele trei traiectorii și interpretați efectul vântului.

%% Problema 1.5 - Traiectoria 3D a unei drone în prezența vântului
clear; clc; close all

R = 12;            % raza manevrei [m]
z0 = 5;            % altitudine inițială [m]
vz = 1.5;          % viteză verticală [m/s]
omega = 0.5;       % pulsație [rad/s]
t = linspace(0,25,500);

vant_x = [0, -2, 2];   % fără vânt, Vest, Est [m/s]
nume = {'fără vânt','vânt spre Vest','vânt spre Est'};

figure
tiledlayout(1,2)

% Vedere 3D
nexttile
hold on
for i = 1:length(vant_x)
    x = R*cos(omega*t) + vant_x(i)*t;
    y = R*sin(omega*t);
    z = z0 + vz*t;
    plot3(x,y,z,'LineWidth',1.5)  % grafic 3D cu functia plot3
end
grid on
xlabel('x - est [m]')
ylabel('y - nord [m]')
zlabel('z - altitudine [m]')
title('Vedere 3D')
view(35,25)
axis tight

% Vedere de sus: proiecția în planul x-y
nexttile
hold on
for i = 1:length(vant_x)
    x = R*cos(omega*t) + vant_x(i)*t;
    y = R*sin(omega*t);
    plot(x,y,'LineWidth',1.5)  % se cere vedere de sus, deci doar in planul xy, functia plot
end
grid on
axis equal
xlabel('x - est [m]')
ylabel('y - nord [m]')
title('Vedere de sus')
legend(nume,'Location','best')

Problema 1.5: graficul 3D, rotit cu mouse-ul3 min

Film în pregătire

Desenarea celor trei traiectorii cu plot3, rotirea graficului și trecerea la vederea de sus, unde se vede deriva produsă de vânt.

Deschide filmul pe pagina lui

Rezultatul rulării

Caz	x(25) [m]	y(25) [m]	z(25) [m]
Fără vânt	11.97	-0.80	42.50
Vânt spre Vest	-38.03	-0.80	42.50
Vânt spre Est	61.97	-0.80	42.50

La t = 25 s, altitudinea este aceeași în toate cele trei cazuri, deoarece modelul consideră vânt numai pe direcția x. Față de cazul fără vânt, poziția finală este deplasată cu 50 m spre Vest, respectiv cu 50 m spre Est.

Cele trei traiectorii ale dronei: vedere 3D și vedere de sus

Interpretare. Fără vânt, proiecția traiectoriei în planul x-y este un cerc, iar în spațiu drona descrie o elice ascendentă. Vântul constant pe direcția x introduce o derivă progresivă: în vederea de sus, cercul nominal este transformat într-o traiectorie de tip trohoidă deplasată spre Vest sau spre Est. Mișcarea pe axa y și viteza de urcare rămân neschimbate în acest model.

PRO 1.6 – Hartă termică a unui modul electronic

Un modul electronic este monitorizat printr-o matrice de 10 x 10 senzori de temperatură. Generați în Matlab o matrice T de 10 x 10 cu valori între 25 și 85 °C (folosiți funcția randi). Afișați harta termică folosind funcțiile imagesc și colorbar. Determinați câte puncte din cele 100 depășesc 70 °C și extrageți pozițiile lor în matrice.

%% PRO 1.6 - Hartă termică (heat map) a unui modul electronic
clear; clc; close all

rng(7)
T = randi([25,85],10,10);

figure
imagesc(T)
colorbar
axis equal tight
xlabel('Coloana senzorului')
ylabel('Linia senzorului')
title('Hartă termică a modulului electronic')

prag = 70;
zone_fierbinti = T > prag;
numar_puncte = sum(zone_fierbinti,'all');
[linie, coloana] = find(zone_fierbinti);

fprintf('Număr puncte peste prag: %d\n', numar_puncte)
TabelZone = table(linie, coloana, T(zone_fierbinti), ...
    'VariableNames', {'Linie','Coloana','Temperatura'});
disp(TabelZone)

Harta termică a matricei de 10 x 10 senzori, cu scala de culori

Rezultatul rulării

Număr puncte peste prag: 23

Linie	Coloană	Temperatura [°C]
1	7	83
1	8	76
1	10	73
2	1	72
2	2	74
2	3	82
2	10	82
4	9	72
4	10	78
5	1	84
5	3	80
5	9	73
6	2	80
6	4	76
7	4	71
7	7	83
7	10	73
8	5	80
8	7	82
8	9	72
9	2	81
9	7	76
9	9	79

Interpretare. 23 senzori depășesc pragul de 70 °C. Tabelul indică poziția fiecărui punct fierbinte (peste valoarea de prag indicată) și temperatura corespunzătoare.

PRO 1.7 – Matrice de senzori pentru o platformă robotică

O platformă robotică are o matrice de senzori S \in \mathbb{R}^{8 \times 8}, care măsoară o mărime normalizată între 0 și 1. Să se determine valoarea maximă, valoarea minimă și poziția senzorului cu valoarea maximă. Dacă maximul apare de mai multe ori, se vor afișa toate pozițiile.

%% PRO 1.7 - Matrice de senzori pentru o platformă robotică
clear; clc; close all

rng(4)
S = rand(8,8);

valoare_maxima = max(S,[],'all');
valoare_minima = min(S,[],'all');

[linie_max, coloana_max] = find(S == valoare_maxima);

fprintf('Valoare maximă: %.4f\n', valoare_maxima)
fprintf('Valoare minimă: %.4f\n', valoare_minima)

table(linie_max, coloana_max, ...
    'VariableNames', {'Linie','Coloana'})

Rezultatul rulării

Valoare maximă: 0.9834
Valoare minimă: 0.0062

Linie	Coloană
6	2

Interpretare. Pentru matricea generată cu rng(4), valoarea maximă este 0.9834 și apare la poziția (6, 2), iar valoarea minimă este 0.0062. Utilizarea funcției find permite identificarea tuturor pozițiilor maximului dacă aceeași valoare apare de mai multe ori.

Întrebări de verificare

Care este diferența dintre ^ și .^ în MATLAB?
De ce este necesară conversia vitezei din km/h în m/s în Problema 1.1?
Ce tip de variabilă este conditie în Problema 1.2?
Care este diferența dintre operatorii & și |?
Ce returnează funcția find atunci când este aplicată unui vector logic?
De ce varianta cu indexare logică din Problema 1.4 este mai eficientă decât varianta cu for?
Ce rol are comanda hold on în reprezentările grafice?
Cum se poate modifica Problema 1.5 pentru a include și o componentă a vântului pe axa y?
Ce scop are utilizarea funcției find la problema 1.5?
Cum se poate salva automat o figură MATLAB într-un fișier PNG?

Temă pentru acasă

Pornind de la Problema 1.5, modificați scriptul astfel încât utilizatorul să poată introduce de la tastatură mărimile care cuantifică acțiunea vântului pe axele x și y. Reprezentați pe același grafic 3D traiectoria nominală și traiectoria afectată de vânt și calculați distanța orizontală dintre pozițiile finale.