Matrisberäkningar med för STM32?

[Rick81]

Nu dribblade Al bort mig igen.....

Var inte hela poängen att köa matris ekvationer på samma enhet som gör parallel 3 ADsamplingar?

Om du tycker koden blir otydlig kan du ju bara använda olika pekare med korrekt namn som pekar på samma area? Pekare tar knappt nåt RAM. Har väldigt svårt att se att du behöver mer än 12 st 144x144 samtidigt i minnet. Och om du bara gör ett eller tre sampel per sekund kommer det ju ta väldigt lång tid innan du kan köra din reglering?

Och hur stor skillnad gör 36x36, 72x73 mot 144x144 matriser i praktiken?

Men sen kan det ju vara så att jag missförstått hela projektet.....

[bearing]

Jag tror/gissar att de flesta cellerna i dom här matriserna har med "machine learning" att göra.
Några få celler fylls med saker som insignal, integral av insignal, derivata, andraderivata, o.s.v. Och samma med utsignalen. Resten av cellerna innehåller parametrar, estimeringar, prediktioner, o.s.v. Eller hur skulle du förklara det, Al?

Det jag inte kan förstå är hur systemet kunde växa från 36x36 till 144x144. Det är alltså 16 gånger mer data? för samma insignal och utsignal?

Inom bilindustrin brukar varje del i mjukvaran få en budget vad gäller minne och processorkraft. Som brukar vara en oerhört stram budget, vilket betyder att ingenjörerna behöver vara kreativa och försöka få tillräckligt bra resultat med en förenklad model av systemet. Man jobbar i matlab, och jämför en "fullständig modell" med olika steg av förenklingar, tills man hittar några bra kompromisser. Sen kollar man hur mycket resurser de olika kompromisserna tar, och om det går att optimera koden på något vis. Jag tycker att det är lite bakvänt att planera för en enorm processor bara för att man någon gång i matlab testar med en stor matris.

[Al_Bundy]

Bygga eget vinner man inget på. Om inte raspbian/ubuntu eller realtimepi funkar så är det bara "bare metal" som gäller.

Tack! Du sparade många timmar åt mig Då blir det Debian server. Tryggt och säker.

[Al_Bundy]

Nu dribblade Al bort mig igen.....

Var inte hela poängen att köa matris ekvationer på samma enhet som gör parallel 3 ADsamplingar?

Om du tycker koden blir otydlig kan du ju bara använda olika pekare med korrekt namn som pekar på samma area? Pekare tar knappt nåt RAM. Har väldigt svårt att se att du behöver mer än 12 st 144x144 samtidigt i minnet. Och om du bara gör ett eller tre sampel per sekund kommer det ju ta väldigt lång tid innan du kan köra din reglering?

Och hur stor skillnad gör 36x36, 72x73 mot 144x144 matriser i praktiken?

Men sen kan det ju vara så att jag missförstått hela projektet.....

Alltså jag håller på att göra machine learning, dvs adaptiv reglering av ett dynamiskt system, Då kan det vara ett behov av matriser.

Men man kan göra det på enklare sätt också. Men boven i dramat stavas Singular Value Decomposition. Det är det bästa matrisverktyget någonsin, men den kostar stora matriser.

Jag tror/gissar att de flesta cellerna i dom här matriserna har med "machine learning" att göra.
Några få celler fylls med saker som insignal, integral av insignal, derivata, andraderivata, o.s.v. Och samma med utsignalen. Resten av cellerna innehåller parametrar, estimeringar, prediktioner, o.s.v. Eller hur skulle du förklara det, Al?

Endast insignal och utsignal mäts för att estimera en matematisk modell. Sedan så används modellen för att prediktivt räkna ut de bästa signalerna.

Det jag inte kan förstå är hur systemet kunde växa från 36x36 till 144x144. Det är alltså 16 gånger mer data? för samma insignal och utsignal?

För jag krävde mer data.

Inom bilindustrin brukar varje del i mjukvaran få en budget vad gäller minne och processorkraft. Som brukar vara en oerhört stram budget, vilket betyder att ingenjörerna behöver vara kreativa och försöka få tillräckligt bra resultat med en förenklad model av systemet. Man jobbar i matlab, och jämför en "fullständig modell" med olika steg av förenklingar, tills man hittar några bra kompromisser. Sen kollar man hur mycket resurser de olika kompromisserna tar, och om det går att optimera koden på något vis. Jag tycker att det är lite bakvänt att planera för en enorm processor bara för att man någon gång i matlab testar med en stor matris.

Med tanke på att industrin använder reglerteknik från 80-talet, medan jag tänker använda reglerteknik från 90-talet.

[Al_Bundy]

Det är ingen annan som håller på med systemidentifiering? I så fall så kan ni söka upp ordet "System identification". Det är ett extremt teoretiskt ämne som jag försöker få till praktik.

Grejen med den moderna systemidentifiering är att den använder sig av Singular Value Decomposition.

\(A = U\Sigma V^T\)

Med vanlig systemidentifiering så löser man system mer exakt t.ex.
\(Ax=B\)

Där \(A \in \Re^{nxn}\)

Därför kan man lösa systemet igenom inversen
\(x = A^{-1}B\)

Men om \(A \in \Re^{nxm}\) så måste man ta pesudoinvers
\(x = (A^T A)^{-1}A^TB\)

Men den nya epoken av systemidentifiering använder sig av SVD
\(x = V \Sigma ^{-1}U^TB\)

Varför? Jo, för den är säkrare och stabilare, men den kräver mer matrisberäkningar.

[bearing]

Mycket beräkningar som görs med matriser blir redundanta om de körs i varje cykel. Man kan ofta förenkla genom att t.ex endast beräkna några diskreta integraler av önskade värden under varje cykel. Och sedan när man integrerat under lämplig tid, sätter man in värdena i sin matris, beräknar inversen, och får en lösning.

Det du gör verkar ju mycket mer avancerat, och du behöver en ett resultat (utsignal) varje cykel, så detta är kanske inte applicerbart, men det bör ju finnas med i bakhuvudet att det kan vara möjligt.

Man använder definitivt matrismatematik i bilar. Men adaptiva algoritmer är inte helt accepterat, då det är mycket viktigt med robusta system, som är verifierbara. Hur ska man kunna verifiera att en adaptiv algoritm alltid anpassar sig på rätt sätt, även efter 10 års drift?

[Al_Bundy]

Mycket beräkningar som görs med matriser blir redundanta om de körs i varje cykel. Man kan ofta förenkla genom att t.ex endast beräkna några diskreta integraler av önskade värden under varje cykel. Och sedan när man integrerat under lämplig tid, sätter man in värdena i sin matris, beräknar inversen, och får en lösning.

Det du gör verkar ju mycket mer avancerat, och du behöver en ett resultat (utsignal) varje cykel, så detta är kanske inte applicerbart, men det bör ju finnas med i bakhuvudet att det kan vara möjligt.

Man använder definitivt matrismatematik i bilar. Men adaptiva algoritmer är inte helt accepterat, då det är mycket viktigt med robusta system, som är verifierbara. Hur ska man kunna verifiera att en adaptiv algoritm alltid anpassar sig på rätt sätt, även efter 10 års drift?

Därför används SVD, eller QR. Tidigare så har adaptiva regulatorer endast passa ett typ av beteende hos ett dynamiskt system t.ex. stokastiskt. Men den moderna identifieringen är dynamisk och kräver inget speciellt beteende just på grund utav att den använder SVD eller QR.

Robusta system är väldigt säkra, men dessa kostar mycket pengar att just få robust. Om du tänker på ämnet Robust Control? Inte mitt favoritämne då jag ämnet är sjukt teoretiskt och skärmad från verkligheten.

Man kan ju ställa sig frågan "Om adaptiva regulatorer är så bra, varför finns dem inte överallt då?". För att dem kostar för mycket. Man måste ha typ en PC för att det ska fungera. Hur länge en dator håller, vet ingen. Men desto mindre saker i en dator, desto mindre kan krångla också.

[Icecap]

Nja - en modern 32-bitars µC motsvarar drygt en 80386!
Den jag använder är "bara" en 50MHz och efter vad jag vet gjordes 386'an inte till högre klockfrekvens än 33MHz.

Så nej, det behövs inte en PC/"dator", det behövs bara att man vet vad man behöver av RAM och PROM samt att man kan få en idé om vilken minsta hastighet som behövs för att klara målen.

Detta är anledningen till att du har blivit frågad om hur många samtidiga matriser du rent faktisk behöver - men dina svar har ju varit ganska dimmiga om detta, sannolikt för att du inte har överblick.

Det är inget problem att få tag i en 32-bitars µC med mer än minne nog och inbyggd FPU samt redigt med hastighet - men det kan kräva en del att starta upp från grunden, mycket datablad och tid för att starta de olika saker som måste köra i botten, typ UART för kommunikation, A/D-omvandlare för datasampling, Timers för de saker som behöver det osv.

[Klas-Kenny]

"Nja - en modern 32-bitars µC motsvarar drygt en 80386!"

En 80386 har väl inte mycket att sätta emot en modern µC..
80386@33MHz fixar ca 11,4 MIPS enligt Wikipedia

Medan på tex. Cortex M4 har man 1,25 DMIPS/MHz, så en sådan på 120MHz (som inte ens är vidare högt för en modern µC) ger 150 DMIPS, alltså mer än tio gånger så snabb som en 80386.


</OT>

[guckrum]

Kod: Markera allt

Robusta system är väldigt säkra, men dessa kostar mycket pengar att just få robust. Om du tänker på ämnet Robust Control? Inte mitt favoritämne då jag ämnet är sjukt teoretiskt och skärmad från verkligheten. 

Du tycker att design av robusta system är skärmat från verkligheten?

I praktiken, mig veterligen, undviker man ofta generella adaptiva algoritmer
för att lösa specifika problem. Åtminstone så länge man inte har förstått
problemet och algoritmen ordentligt och kan säkerställa att lösningen
uppfyller de krav man har. (Krav kan också vara på saker som
observerbarhet och förståelse.) Generella lösningar kan vara en bra grogrund
till en lösning och kan användas till att jämföra sig mot, till exempel. Och,
som sagt, ibland är de också en utmärkt lösning.

Att använda så kallade "black box"-lösningar innebär bland annat att man inte
kan göra ordentlig riskanalys och att man inte kan säkra beteende eller kanske
ens konvergens. I ditt systemidentifieringsfall kommer du att sprida den entropi
du har över ett stort antal matriselement, och frågan är vilka av dem som är
signifikanta. Adderar du rekursion till detta kan du få alla möjliga intressanta
effekter på grund av saker som numerik och rank, vilket kan leda till överraskande
placeringar av systemets poler, till exempel. I vissa fall går det utmärkt, i andra
mindre bra. Frågan är vilka fall som fungerar.

Att förstå problemet/systemet och därefter välja en specifik lösning är naturligtvis
jobbigare och kräver mera kunskap och erfarenhet, men risken att bli biten
kommer att vara mycket mindre.

[Al_Bundy]

Jag håller fortfarande på att utveckla min MPC nu.

Så om ni väntar ett tag så ska jag ge ett ytterligare exempel program.

Då kan vi se exakt hur många matriser jag använder

[Al_Bundy]

Kod: Markera allt

Robusta system är väldigt säkra, men dessa kostar mycket pengar att just få robust. Om du tänker på ämnet Robust Control? Inte mitt favoritämne då jag ämnet är sjukt teoretiskt och skärmad från verkligheten. 

Du tycker att design av robusta system är skärmat från verkligheten?

I praktiken, mig veterligen, undviker man ofta generella adaptiva algoritmer
för att lösa specifika problem. Åtminstone så länge man inte har förstått
problemet och algoritmen ordentligt och kan säkerställa att lösningen
uppfyller de krav man har. (Krav kan också vara på saker som
observerbarhet och förståelse.) Generella lösningar kan vara en bra grogrund
till en lösning och kan användas till att jämföra sig mot, till exempel. Och,
som sagt, ibland är de också en utmärkt lösning.

Att använda så kallade "black box"-lösningar innebär bland annat att man inte
kan göra ordentlig riskanalys och att man inte kan säkra beteende eller kanske
ens konvergens. I ditt systemidentifieringsfall kommer du att sprida den entropi
du har över ett stort antal matriselement, och frågan är vilka av dem som är
signifikanta. Adderar du rekursion till detta kan du få alla möjliga intressanta
effekter på grund av saker som numerik och rank, vilket kan leda till överraskande
placeringar av systemets poler, till exempel. I vissa fall går det utmärkt, i andra
mindre bra. Frågan är vilka fall som fungerar.

Att förstå problemet/systemet och därefter välja en specifik lösning är naturligtvis
jobbigare och kräver mera kunskap och erfarenhet, men risken att bli biten
kommer att vara mycket mindre.

Det är förvånade att många piskar på mig att jag inte läser det andra skriver, men hade du läst mitt citat du citerade så förstår du vad jag menar.

Robusta system är väldigt säkra, men dessa kostar mycket pengar att just få robust. Om du tänker på ämnet Robust Control? Inte mitt favoritämne då jag ämnet är sjukt teoretiskt och skärmad från verkligheten.

De flesta som håller på med design av reglersystem använder sig i praktiken alltid av PID regulatorer. Det har med att PID fungerar och PID gör sina jobb. Men sedan finns det ett ämne, som jag talade om i mitt citat, som heter Robust Control.

Robust Control handlar om att forma kretsförstärkningen. Jag har läst lite om Robust Control och det är enormt skärmad från verkligheten då fokus ligger endast på teori, utan att ta hänsyn till verklighet, vilket Robust Control ska egentligen vara. Om vi ska bygga en robot här på Jorden och den ska till Mars om 6 månader, då kanske vi ska teoretiskt designa en styrning och balans för roboten, grundat på generella fysiska formler. Här tar vi inte hänsyn till vinden som blåser, eller gravitation. Vi försöker bara vara så robust som möjligt igenom att bygga upp en generell modell och teoretiskt vrida och vända på kretsförstärkningen.

I praktiken används Robust Control inte. Den är för ny och obeprövad. Dessutom kan en P-regulator eller en PID bli minst lika robust. Så ur en akademisk synpunkt så är Robust Control intressant, men ur en verklighets vinstdrivande synpunkt så handlar det om att använda en lämplig metod.

[Al_Bundy]

Nu har jag skrivit två program. Ena programmet estimerar en matematisk modell av mätdata och det andra programmet räknar ut vilka insignaler man ska ha för att uppnå ett börvärde, dvs en referens.

Här kan ni se hur många matriser jag använder mig utav. Det är alltså många. Nu fattas det bara att jag kopplar ihop dessa två och sedan låter jag utsignalen U styra vilken PWM frekvens det ska vara.

Se bifogad fil. Ny uppdatering. Denna gång är ALLT double!

Identifieringsprogrammet:

Kod: Markera allt

#include <time.h>
#include "LinearAlgebra/declareFunctions.h"


int main() {

	/*
	 * G(s) = 1/(s^2 + 1s + 3)  - Model
	 * y = measured output values
	 * u = measured input values
	 */

	clock_t start, end;
	float cpu_time_used;
	start = clock();

	double y[144] = { 0.00000, 0.49525, 1.43863, 2.13779, 2.30516, 2.05713,
			1.69220, 1.45608, 1.42777, 1.54146, 1.67927, 1.75624, 1.75400,
			1.70478, 1.65394, 1.62996, 1.63549, 1.65594, 1.67426, 1.68125,
			1.67752, 1.66930, 1.66285, 1.66102, 1.66300, 1.66621, 1.66842,
			1.66880, 1.66786, 1.66664, 1.66591, 1.66588, 1.66629, 1.66675,
			1.66698, 1.66695, 1.66678, 1.66661, 1.66654, 1.66657, 1.66664,
			1.66670, 1.66672, 1.66670, 1.66667, 1.66665, 1.66665, 1.66666,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66666, 1.66666,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667, 1.66667,
			1.66667, 1.66667, 1.66667, 1.66667, 1.66667 };

	double u[144] = { 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
			5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
			5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
			5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
			5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
			5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
			5, 5, 5, 5, 5, 5, 5, 5, 5 };

	// Toeplitz matrix
	double toe[144 * 144];
	toeplitz(u, toe, 144);

	// Upper triangular
	double tru[144 * 144];
	triu(toe, tru, 0, 144, 144);

	// inverse
	inv(tru, 144);

	// Multiplication
	double g[144];
	mul(y, tru, false, g, 1, 144, 144);

	// Create hankel
	double H0[144 * 144];
	hankel(g, H0, 144, 1);
	double H1[144 * 144];
	hankel(g, H1, 144, 2);

	// Cut hankel into half
	double H0_half[72 * 72];
	double H1_half[72 * 72];
	cut(H0, 144, 144, H0_half, 0, 71, 0, 71);
	cut(H1, 144, 144, H1_half, 0, 71, 0, 71);

	// Do SVD
	double U[72 * 72];
	double S[72 * 72];
	double V[72 * 72];
	svd(H0_half, U, S, V, 72, 72);

	// Model reduction to second order model
	double Un[72 * 2];
	double Sn[2 * 2];
	double Vn[72 * 2];
	cut(U, 72, 72, Un, 0, 71, 0, 1);
	cut(S, 72, 72, Sn, 0, 1, 0, 1);
	cut(V, 72, 72, Vn, 0, 71, 0, 1);

	// Create A, B, C
	double Sn_sqrt_negative[2 * 2];
	double Sn_sqrt_positive[2 * 2];
	copy(Sn, Sn_sqrt_negative, 2,2);
	copy(Sn, Sn_sqrt_positive, 2,2);
	diagpower(Sn_sqrt_negative, -0.5, 2, 2);
	diagpower(Sn_sqrt_positive, 0.5, 2, 2);

	/*
	 * C = Un*Sn^(1/2);
	 * Cd = C(1, 1:2)
	 */
	double C[72 * 2];
	double Cd[1 * 2];
	mul(Un, Sn_sqrt_positive, false, C, 72, 2, 2);
	cut(C, 72, 2, Cd, 0, 0, 0, 1);

	/*
	 * Ad = Sn^(-1/2)*Un'*H1*Vn*Sn^(-1/2);
	 */
	double A[72 * 2];
	double A1[72 * 2];
	double Ad[2 * 2];
	double Ad1[2 * 2];
	mul(Vn, Sn_sqrt_negative, false, A, 72, 2, 2);
	mul(H1_half, A, false, A1, 72, 72, 2);
	tran(Un, 72, 2);
	mul(Un, A1, false, Ad1, 2, 72, 2);
	mul(Sn_sqrt_negative, Ad1, false, Ad, 2, 2, 2);

	/*
	 * B = Sn^(1/2)*Vn'
	 * Bd = B(:, 1);
	 */
	double B[2 * 72];
	double Bd[2 * 1];
	tran(Vn, 72, 2);
	mul(Sn_sqrt_positive, Vn, false, B, 2, 2, 72);
	cut(B, 2, 72, Bd, 0, 1, 0, 0);

	/*
	 * Print A, B, C
	 */
	printf("A Matrix: \n");
	print(Ad, 2, 2);

	printf("B Matrix: \n");
	print(Bd, 2, 1);

	printf("C Matrix: \n");
	print(Cd, 1, 2);

	end = clock();
	cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
	printf("\nTotal speed  was %f,", cpu_time_used);
	return 0;
}

Kontrolprogrammet

Kod: Markera allt

#include <time.h>
#include "LinearAlgebra/declareFunctions.h"

/*
 * Here you can set the sizes for the matrices
 */

#define row_a 2 // A row
#define column_a 2 // A column

#define row_b 2 // B row, the same row as A.
#define column_b 1 // B column, the same column as D

#define row_c 1 // C row, the same row as D
#define column_c 2 // C column, the same column as A

/*
 * Create a state space model
 */
double A[row_a * column_a] = { 0, 1,
							  -3, -2};

double B[row_b * column_b] = { 0,
								1
};

double C[row_c * column_c] = {1, 0};

double D[row_c * column_b] = {0};

double X[row_a] = {0, 0};

double R =  6;

// Create the length of the observability matrix, notice it will have the dimension (row_o * row_c + row_c) x column_b
#define row_o 30

// Create the widt of the lower triangular toeplitz H matrix, notice that it will have the dimension (row_o * row_c + row_c) x (column_h * column_b)
#define column_h 29 // column_h < row_o - always!

// Define the column of the reference vector - Standard is 1
#define column_ry 1

int main() {

	/*
	 * Model predictive control
	 */

	clock_t start, end;
	float cpu_time_used;
	start = clock();

	/*
	 * Create the Observabillity matrix and the
	 */

	double O[(row_o * row_c) * row_a];
	double O_[(row_o * row_c) * row_a]; // This is for the lower triangular toeplitz matrix
	double A_[row_a * column_a];
	double C_[row_c * column_c];

	for (int i = 1; i <= row_o; i++) {
		copy(A, A_, row_a, column_a); // Copy A -> A_
		mpower(A_, row_a, i); // Power A_^i
		mul(C, A_, false, C_, row_c, column_c, column_a); // C_ = C*A_
		insert(C_, O, row_c, column_c, row_a, (i-1)*row_c, 0); // Insert O = [CA^1; CA^2; CA^3; ... ; CA^row_o];

		copy(A, A_, row_a, column_a); // Copy A -> A_
		mpower(A_, row_a, i - 1); // Power A_^(i-1)
		mul(C, A_, false, C_, row_c, column_c, column_a); // C_ = C*A_
		insert(C_, O_, row_c, column_c, row_a, (i-1)*row_c, 0); // Insert O_ = [CA^0; CA^1; CA^2; ... ; CA^(row_o-1)];
	}

	//print(O, row_o * row_c, row_a);

	/*
	 * Create the lower triangular toeplitz matrix
	 */
	double H[(row_o * row_c) * (column_h * column_b)];
    zeros(H, row_o * row_c, column_h * column_b);

	// T = O_*B - Observabillity matrix times B
    double T[(row_o * row_c) * column_b];
	mul(O_, B, false, T, row_o * row_c, row_a, column_b);


	// Lower triangular toeplitz matrix of T = [C*A^0*B; C*A^1*B; C*A^2*B; C*A^3*B; ...; C*A^(row_o-1)*B]
	for (int j = 0; j < column_h; j++){
		insert(T, H, row_o * row_c, column_b, column_h * column_b, 0, j*column_b);
		move(T, row_o * row_c, column_b, row_c , 0);
	}

	//print(H, row_o * row_c, column_h * column_b); // H matrix


	/*
	 * Compute U = pinv(H)*(Ry*R - O*X), where R is our reference vector, X is our initial state vector
	 */

	pinv(H, row_o * row_c, column_h * column_b); // Pseudo inverse of H. Using the SVD method

	double Ry[(row_o * row_c)*column_ry];
	ones(Ry, row_o * row_c, column_ry);
	scale(Ry, R/2, row_o * row_c, column_ry); // Ry*R = Ry (R need to be divided with 2)

	double OX[(row_o * row_c)*column_ry];
	mul(O, X, false, OX, row_o * row_c, row_a, column_ry); // O*X

	double Ry_OX[(row_o * row_c)*column_ry];
	sub(Ry, OX, Ry_OX, row_o * row_c, column_ry, column_ry); // Ry-O*X

	double U[(column_h * column_b)*column_ry];
	mul(H, Ry_OX, false, U, column_h * column_b, row_o * row_c, column_ry); // U = pinv(H)*(Ry-O*X);

	/*
	 * Our best input values
	 */
	print(U, column_h * column_b, column_ry);



	end = clock();
	cpu_time_used = ((float) (end - start)) / CLOCKS_PER_SEC;
	printf("\nTotal speed  was %f,", cpu_time_used);
	return 0;
}

[Al_Bundy]

Nu är jag sugen på att börja implementera detta i verkligheten. Jag har ju som sagt skapat två program som ska göra mina beräkningar. Jag funderar på att implementera detta på en enkortsdator. Vad för enkortsdator, spelar ingen roll. Troligtvis blir det en Paj med tanke på priset.

Jag tänkte först att ha mina program på en vanlig PC, men då kom jag på att jag måste ha kommunikation mellan USB och ett t.ex STM32 kort där STM32 kortet fungerar som en I/O modul. Men en Paj har ju redan I/O gångar så jag använder då en Paj. Dessutom så håller väll en Paj riktigt länge, med tanke på att så länge man inte skriver till SD kortet så är väll en Paj odödlig?

Då är min fråga. Hur fungerar RealTimePi?
https://github.com/guysoft/RealtimePi

Är det bara testbaserat?

[guckrum]

Jag hänger nog inte med riktigt, skall y och u vara ut- respektive insignal
till systemet G(s) = 1/(s^2 + 1s + 3)? I så fall borde väl y gå mot noll
stationärt, eftersom G(s) har komplexa konjugatpoler i vänster halvplan?
Eller jag kanske missar något. Vad säger ditt program?

G(s) - om man tänker i Laplacetermer så är det bra att minnas att man
har att göra med ett samplat system som bara "approximerar" ett tids-
kontinuerligt dito. Detta är kanske bara något man behöver ta hänsyn till
vid konvertering mellan domänerna, men bra att ha i bakhuvudet när
man jobbar, hur som helst.

Och robustheten. Som jag förstår detta är det antalet sampel du tar
hänsyn till som styr ordningen på matriserna, och dessa sampel plockas
ut med ett "rektangelfönster". Jag vet inte hur teorin för det här ser ut,
men beroende på när du samplar och kombinationen av insignal, system
och utsignal så kommer approximationen att bli olika "bra"? Och vad
händer egentligen om det adderas brus på signalen? För att vara mera
precis, hur ser fördelningen av felet ut dels för olika ordningar på
matriserna, dels för olika nivåer av, exempelvis, additivt gaussiskt brus?

För mig är det väldigt ointuitivt att "systemidentifiera" ett andra ordningens
system med matriser av ordning 144, så jag är lite nyfiken.